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THEORY  AND  SIMULATION  OF  A  GATED  FIELD  EMITTER: 
ANALYSIS  AND  INCORPORATION 
INTO  MACROSCOPIC  DEVICE  CHARACTERIZATION 


INTRODUCTION 

Recent  improvements  in  gated  field  emission  array  (FEA)  technology  provide  a  new  alternative  to 
thermionic  electron  beam  sources  [1]  for  emission  gating  at  frequencies  above  ultra  high  frequency 
(UHF).  This  new  opportunity,  however,  is  subject  to  the  integration  of  FEA  technology  into  the 
vacuum  tube  environment.  The  principal  advantages  of  microfabricated  FEAs  are  the  negligible  tip- 
to-gate  transit  time  and  the  high  transconductance.  FEA  structures  can  modulate  the  beam  density  at 
high  frequency  and  with  good  spatial  localization,  thereby  making  possible  significant  improvements 
in  amplifier  performance  [2].  They  have  also  been  successfully  used  in  the  creation  of  flat  panel 
displays  [3]  and  other  devices.  Electrons  from  the  ultra  sharp  tips  (radius  of  curvature  is  on  the  order 
of  100  A)  can  be  emission  modulated  by  applying  an  oscillating  voltage  to  the  gate  that  is  coplanar 
with  the  emitter  tips.  FEAs,  therefore,  promise  orders  of  magnitude  improvement  in  frequency 
response  over  the  thermionic  cathodes  (for  which  the  gate  plane  is  on  the  order  of  20  pm  from  the 
cathode  surface)  due  to  limitations  associated  with  transit  time  effects. 

Because  of  this,  considerable  interest  in  developing  an  FEA  cathode  for  inclusion  into  a  twystrode 
or  klystrode  amplifier  configuration  has  developed.  Assumptions  regarding  the  relation  between  the 
cathode  current  and  the  gate  voltage  are  required  to  estimate  inductive  output  amplifier  (lOA) 
performance  characteristics  (e.g.,  gain,  bandwidth,  frequency  of  operation,  efficiency)  [4]. 
Consequently,  the  prediction  of  device  performance  has  been  impeded  by  the  lack  of  a  simple 
methodology  for  comparison  of  disparate  FEA  structures.  Different  fabrication  methods  produce 
widely  different  gate  apertures,  emitter  heights,  radii  of  curvature  at  the  tip,  and  other  geometric,  as 
well  as  material,  differences.  The  FEA  unit  cell,  therefore,  exhibits  different  current-voltage 
characteristics,  transconductances,  and  capacitances. 

Experimental  data  are  well  characterized  by  the  Fowler-Nordheim  (FN)  equation  [5,6],  and 
modeling  efforts  typically  take  advantage  of  this.  However,  numerous  physical  effects  mitigate 
against  simple  first-principles  analysis  of  field  emission  from  an  FEA  tip,  and  therefore  invite  the 
question  “Is  FEA  performance  predictable?”  Effects  such  as  the  influence  of  the  anode  voltage 
during  low  current  operation,  or  space  charge  effects  during  high  current  operation,  may  be 
amenable  to  a  first-principles  analysis,  and  a  fundamental  understanding  of  these  effects  seems 
possible.  Other  effects,  well  known  to  fabricators,  include  contamination,  outgassing,  seasoning, 
migration,  dulling,  and  snapback.  These  are  influenced  by  manufacturing,  conditioning,  and 
operational  history,  and  are  difficult  to  interpret  in  terms  of  first  principles.  Typically,  these  effects 
are  controlled  or  held  constant  by  holding  to  a  strict  operational  regimen.  Another  effect  that  can 
prevent  first-principles  analysis  of  emission  from  an  FEA  tip  is  geometric  irregularity  of  the  tips. 
Recent  transmission  electron  microscope  (TEM)  photographs  indicate  that  actual  tip  shapes  are 
highly  irregular,  whereas  convenient  analytical  models  assume  spherical  or  elliptical  shapes.  Certain 
data  support  the  notion  of  emission  from  a  single,  or  perhaps  a  few,  sites  on  a  tip  [7],  in  marked 
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contrast  to  the  idea  of  emission  occurring  continuously  over  a  broad  area.  Although  the  actual 
quantum-mechanical  process  may  have  a  classical  counterpart,  key  ingredients,  such  as  the  electronic 
charge,  remain  quantized.  At  scales  of  interest,  only  a  single  electron  may  exist  in  the  vicinity  of  the 
gate,  which  may  alter  our  understanding  of  “space  charge”  in  depressing  surface  fields.  “Turn 
on,”  in  which  a  seemingly  dysfunctional  tip  under  some  subtly  changing  influence  suddenly  begins 
to  emit,  has  been  observed,  again  suggesting  lack  of  continuity.  In  light  of  this,  it  appears  improbable 
that  a  truly  predictive  and  comprehensive  theory  or  simulation  analysis  can  be  found. 

Statistical  Models 

Discrete:  The  Linear  Distribution  in  B 

But  before  such  pessimism  is  accepted,  the  fact  remains  that  the  Fowler-Nordheim 
characterization  of  experimental  data  is  of  great  utility,  even  though  the  fabrication,  processing,  and 
even  measurement,  procedures  vary  widely.  We  must  therefore  ascertain  exactly  what  can  be  inferred 
from  the  experimental  data  through  a  statistical  analysis.  Let  the  individual  tips  in  a  field  emission 
array  (FEA)  individually  obey  a  Fowler-Nordheim  relation  (a  similar  analysis  has  been  performed  by 
Levine  [8]  at  lower  emissivities  for  flat  displays).  The  ith  tip  is  characterized  by  (the  “/n”  subscript 
on  A  and  B  is  suppressed): 


/,<V)  =  A,.F2exp(-B,./V^)  ,  (1) 

where  Vg  is  the  gate  voltage  (below,  the  “g”  subscript  is  suppressed).  A  plot  of  ln(/,/Vg2)  vs  MVg  is 
linear.  Below,  we  refer  to  this  as  “Fowler-Nordheim  (FN)  behavior”  (see  Appendix  A). 

Assume  that  the  number  of  tips  characterized  by  B,-  is  n,-.  Then  the  average  current  per  tip  from 
an  FEA  containing  N  tips  is: 

=  N  /?i  =  N  '  ^1  ’  (2) 

where  Ag  and  Bg  remain  to  be  specified,  and  where  L  /i,  =  N.  Consider  the  case  where  the  B’s  are 
linearly  distributed  according  to  Bi  =  Bg  +  (i-l)  AB  /(Ng-1),  which  we  refer  to  as  the  "Linear 

Distribution."  Further,  assume  that  =  (N  /  Nb)  Ag.  Because  increases  with  tip  radius  (as 

does  Bf„),  this  implies  that  a  linear  distribution  in  B’s  does  not  entail  a  uniform  distribution  in  tip 

radius:  «,•  is  larger  for  the  sharper  tips.  Performing  the  summation  in  Eq.  (2),  </>  is  given  by 


1  -  exp 


1  -  exp 


Ns-\  V  I 


1  ab] 

Ng-l  V  j 


(3) 


As  Nb^  °°  (the  continuum  limit),  </>  =  (V/AB)  IgiV)  [l-exp(AB/V)].  The  continuum  limit  may  also 
be  obtained  directly  by  introducing  the  ("linear")  distribution  function  {n^Ai!  NAg  )  =^f(B)  dB  =  dB 
6(B-Bo)  G(Bg  +  AB  -  B)/AB,  and  transforming  the  summation  in  Eq.  (2)  into  an  integral.  f(B)  is  so 
defined  that  its  integral  over  all  B  is  unity. 
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For  A5/V  «1,  we  see  that  </>  =  IgiV),  i.e.,  the  tips  are  highly  uniform.  For  AB/V  sufficiently 
large  (but  small  compared  to  Nb),  </>  =  (V  /  AB)  Ig(V).  Strictly  speaking,  the  coefficient  (V  /  AB) 
will  cause  nonlinearities  in  an  FN  plot,  but  the  question  remains,  to  what  extent  are  these  non- 
linearities  apparent?  We  return  to  this  shortly. 

Continuum 


Consider  the  continuum  extension  of  Eq.  (2)  for  which  the  largest  and  smallest  B,’s  are  denoted 
by  B+  and  B_,  respectively.  Let  the  distribution  of  B’s  be  analogous  to  the  Gaussian  distribution: 


/(5)  = 


exp[-(B-B^)^/(y2] 


exp[-(B-B„)2/a2]r/B 


(4) 


Because  (B+  -  Bg)  need  not  equal  {Bg  -  B_),  Bg  will  not  in  general  correspond  to  the  mean  B  nor  will 
a  correspond  to  the  standard  deviation.  At  present,  they  are  simply  parameters  specifying  the 
distribution  of  B’s.  In  what  follows,  we  do  not  let  B±  go  to  ±  «» immediately.  For  example,  consider 
the  characteristics  of  the  tips  created  by  laser  holography  [12]:  the  radius  of  curvature  is  on  the  order 
of  10  to  50  A.  We  expect,  given  the  size  of  individual  atoms,  that  B_  will  likewise  be  constrained  to 
this  range,  although  B^  need  not  be.  Consequently,  the  distribution  of  B’s  for  these  tips  is  probably 
not  a  full  Gaussian,  and  we  therefore  do  not  approximate  the  denominator  of  Eq.  (4)  by  Vn  a. 
Evaluating  </>  =  Jf(B)  l(V)  dB  gives: 


(/)  =  /^(V)exp(772) 


Erf(x^.  +  T])-  Erf(jt_  +  77) 
Erf(x+)-Erf(x_) 


(5) 


where  =  (B^-Bg)/c  and  rj  =  <jl(2V),  and  where  Erf(A:)  =  (2/Vjt){J(j“  exp(-;c2)  dx]  is  the  error 
function  [11]. 


In  the  limit  that  Xo  =  0  and  A  -  (x+-  x^)/2  approaches  «»,  we  find: 

(/)  =/„(F)exp(Ti2)(i  -AATi2exp(-A2) 


(6) 


which  corresponds  to  the  findings  of  Levine:  FN-like  behavior  will  occur  only  if  a  is  small,  where  a 
is  approximately  the  standard  deviation.  Physically,  there  is  not  a  preponderance  of  tips  with  the 
smallest  value;  therefore,  linearity  in  the  FN  behavior  restricts  the  distribution  of  B  values  that  may 
be  present. 

Second,  consider  the  case  where  A  is  small  by  comparison  to  x:^±T|  and  jc^,  approximating  the 
linear  distribution  of  5’s  encountered  in  Eq.  (3)  (note  that  need  not  equal  BJo),  It  can  be  shown 
that 


(/)  ^IJiy)e-'^o^.  (7) 

The  combination  2jc^ti  =  (5+  +  -  2B^)I2V\  consequently,  <I>lV^  =  exp(-Bf^/V),  where  B^  =  {B+  + 

B_)/2,  similar  to  the  discussion  following  Eq.  (3), 
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Finally,  consider  the  case  where  B_  -  and  B^  =  This  case  resembles  half  of  the  Gaussian 
distribution  (and  is  referred  to  as  the  "half-Gaussian"  distribution).  The  current  becomes 

(/)  =  /^(V)exp(Ti2)[l-Erf(Ti)] 

or  7]  large,  1-Erf (rj)  -  txp(-ri^)/ri^7r  ,  and  </>  -  (2V/a)  Again,  a  coefficient  linear  in  V  is 

obtained. 


Least‘Squares  Fit  of  <I>:  FN  -  Linearity  of  a  Distribution  of  Emitters 


Let  us  define  y  =  ln(</>/V^j  and  x  =  //V  (which  does  not  correspond  to  the  x  of  the  previous 
section).  We  desire  to  find  the  best  polynomial  approximation  to  y,  we  denote  as  y'.  We  define  Xq  = 
(^max  +  ymin)^max^min^  similarly,  S  =  (Vmax  ”  ^min)^max^ min^  where  Vmax  ^^d  Vmin  denote 
the  range  over  which  we  desire  the  best  fit.  The  polynomial  approximation  of  order  n,  which  is  the 
best  least-squares  approximation,  can  then  be  given  by 


(8a) 


Pk(x)  is  the  Legendre  polynomial  of  degree  k,  and  the  Q  are  determined  by 

=  (^+ 1)  J_  j  3'W  dz 


(8b) 


where  we  have  introduced  z  =  (x-Xo)/S.  The  Error  and  Coefficient  of  Correlation  are  analogous  to 
the  equivalent  formulae  in  Appendix  A.  For  a  best  linear  fit,  we  truncate  the  summation  in  Eq.  (8a) 
at «  =  1.  As  shown  in  Appendix  A,  the  presence  of  V-dependent  coefficients  does  alter  the  estimates 
of  the  best  Afn  and  Bj^.  Consequently,  the  nature  of  the  distribution  of  the  emitters  as  a  function  of 
Bfii  (or,  equivalently,  as  a  function  of  tip  radius)  can  greatly  affect  the  estimate  of  the  array’s  Afr  and 
Bfn  parameters,  as  compared  to  those  parameters  for  a  particular  FEA  unit  cell.  In  general,  a 
distribution  of  emitters  will  tend  to  lower  the  estimate  of  Afn  and  raise  the  estimate  of  Bfn  compared  to 
the  unit  cell  by  an  amount  depending  on  the  nature  of  the  V-dependence  of  the  coefficient  of 
The  exact  distribution  of  emitters  for  a  given  array  of  field  emitters  is  not,  in  general,  known.  In  the 
absence  of  tip-by-tip  examination,  it  is  difficult  to  estimate,  thereby  motivating  our  restriction  of  the 
cases  of  interest  to  those  where  an  analytic  solution  is  obtainable,  namely,  the  linear  and  half-Gaussian 
distributions. 


In  the  absence  of  a  tip-by-tip  examination  of  an  FEA  or  scaling  experiments,  and  the  nature  of 
the  approximations  used  in  a  statistical  analysis,  the  linearity  of  experimental  data  appears  to  support 
arguments  that  the  tips  are  very  uniform,  or  that  emission  is  due  to  one  kind  of  tip  or  protrusion.  A 
statistical  analysis  coupled  with  a  comprehensive  modeling  of  the  FEA  unit  cell  offers  some  hope  that 
analysis  of  FEA  tips  from  a  Fowler-Nordheim  approach  may  offer  practical  understanding  of  some 
aspects  of  their  operation,  and  that  geometric  irregularity,  whether  it  exists  or  not,  does  not  prevent 
such  analyses.  Certainly,  the  one-dimensional  Fowler-Nordheim  equation  relating  current  density  to 
field  at  the  surface  (with  modification  to  account  for  three-dimensional  effects)  is  to  a  large  extent 
valid,  and,  as  we  argue  and  show  below,  does  not  preclude  (experimental)  linear  total  current  vs  gate 
voltage  relations  for  an  FEA.  Considerable  light  may  be  shed  on  this  subject  through  electron  beam 
profile  characterizations  of  individual  or  small  numbers  of  tips;  this  is  currently  underway  at  NRL 
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[13].  Regardless  of  whether  high  uniformity  or  small  tip  domination  is  responsible,  the  approach 
developed  here  addresses  either  situation. 

One  final  point  to  make  is  that,  due  to  the  behavior  of  the  current  density  J{F)  with  respect  to 
field,  the  total  current  from  a  tip  will  be  dominated  by  those  regions  in  which  the  field  is  largest. 
Therefore,  geometric  detail  and  variations  will  tend  to  be  obscured  by  the  sites  with  the  sharpest 
curvature.  The  fact  that  the  1{V)  characteristics  are  dominated  by  the  dependence  on  one  parameter 
(V)  suggests  that  a  “peak  field”  is  responsible  for  most  of  the  emission,  i.e.,  the  emission  appears  to 
be  independent  of  the  precise  field  distribution.  This  observation  is  reinforced  in  subsequent 
derivations. 

APPLICATION  OF  1-D  TRANSPORT  EQUATIONS  TO  3-D  GEOMETRIES 

The  most  widely  used  equation  to  model  electron  field  emission  is  the  Fowler-Nordheim 
equation.  This  relates  current  density  J  to  applied  field  F  according  to  the  relation  (Appendix  B) 

Jfn^T)  =  af„{F)  (9) 

where  the  field  dependence  of  ay>,  and  bf„  is  shown  explicitly,  and  is  due  to  the  field  dependence  of 
the  function  v(y)  and  t(y)  (Appendix  B).  (For  convenience,  the  tables  are  grouped  at  the  end  of  the 
text.)  Implicitly  assumed  in  Eq.  (9)  is  the  field  emission  potential  barrier,  with  an  image  charge 
contribution  of  the  form 


V(x)  =  x-Fx-f,  (10) 

where  %  =  (()  +  |i  in  the  case  of  metals,  or  the  electron  affinity  in  the  case  of  semiconductors.  The  use 
of  the  (Q/x)  form  of  the  image  charge  presupposes  emission  from  a  planar  surface.  In  reality,  surface 
curvature  affects  the  form  of  the  image  charge  contribution  [14-16],  and  a  careful  derivation  shows 
that  the  form  of  the  denominator  should  be  (jc  +x:^)  rather  than  merely  (x)  [17,18].  Furthermore,  the 
form  of  Eq.  (9)  requires  modification  for  semiconductors  [19],  although  for  present  purposes  we 
assume  metallic  parameters. 

Its  limitations  notwithstanding,  Eq.  (9)  is  often  applied  to  field  emission  from  multidimensional 
structures.  This  is  done  by  making  the  simple  replacements  I  =  bJ  and  F  =  6 (a  procedure  that  has 
been  roundly  criticized  [15,20],  but  nevertheless  widely  used)  to  “derive”  current  voltage  relations 
amenable  to  experimental  interpretation  of  data  of  the  form 


/(Vp=A^V^e-®/'>''«,  (11) 

where  and  are  assumed  to  be  constant.  Two  difficulties.  First,  the  form  of  Eq.  (9)  does  not 
ensure  that  A^„  and  are  even  constant  with  respect  to  Vg.  Other  curves  are  also  linear,  to  a  good 
approximation,  on  an  FN  plot  (Appendix  A).  Second,  Af„  and  are  experimental  parameters 
obtained  from  the  fitting  of  data;  they  cannot  be  fundamentally  derived  using  Eqs.  (9)  and  (10) 
(especially  in  the  absence  of  a  statistical  analysis)  due  to  the  ad  hoc  nature  of  {b)  and  (B):  the  problem 
of  “measuring”  J  and  F  has  simply  been  deferred  to  “measuring”  b  and  6.  Nevertheless,  the 
comparison  of  FEA  tips  and  geometries  requires  some  figure  of  merit,  and  as  will  be  shown,  the  FN  A 
and  B  parameters  extrapolated  from  the  portion  of  the  /(Vg)  curve  where  the  device  is  intended  to 
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operate  serves  that  purpose.  Consequently,  there  is  a  pressing  need  to  make  some  analytical  or 
simulation-based  estimation  of  and  Bf„.  In  the  past,  they  have  been  estimated  by  using  finite- 
element  methods  or  other  techniques  [20-24],  But  this  is  computationally  intensive  and  does  not 
transparently  reveal  the  parametric  dependence  of,  for  example,  the  I(V)  relations,  on  geometry  and 
materials.  This  report  provides  a  simple  analytical  model  (Saturn)  ,  explores  its  limitations,  compares 
it  with  a  full  numerical  model  based  on  the  boundary  element  method  [25,26],  and  proposes  a  semi- 
numerical  model  that  endeavors  to  overcome  the  limitations  of  the  Saturn  model,  but  nevertheless  is 
computationally  far  less  intensive  than  the  boundary  element  model.  The  analytical  techniques  used 
by  the  Saturn  model  borrow  heavily  from  earlier  analytical  work  dealing  with  field  emission  tips  in  a 
diode  configuration  [27]. 

The  simple,  or  Saturn,  model,  replaces  the  field  emission  tip  with  a  sphere  and  the  gate  with  a 
charged  ring.  An  analytic  relation  then  exists  between  the  gate  radius  and  voltage,  the  sphere  radius, 
the  field,  and  its  angular  variation  on  the  sphere,  from  which  the  total  emitted  current  can  be  derived 
analytically.  For  reasons  to  be  discussed,  the  estimation  of  gate  voltage  is  problematic  in  this  model, 
but  the  utility  of  the  model  lies  more  in  its  ability  to  makes  qualitative  predictions  that  are  borne  out 
in  the  boundary  element  treatment,  rather  than  as  a  precise  computational  tool.  Finally,  the 
methodology  used  in  the  simple  model  motivates  the  approximations  made  in  the  semi-numerical 
approach.  The  boundary  element  method  models  the  tip,  gate,  and  anode  surfaces  with  charged 
ribbons.  At  present,  attention  is  restricted  to  electrodes  with  azimuthal  symmetry.  The  charge  on  each 
ribbon  is  determined  by  a  matrix  solution  of  Poisson’s  equation  for  a  given  tip,  gate,  and  anode 
potential.  The  field  and  area  of  each  ribbon  is  then  calculated,  from  which  the  current  per  ribbon  can 
be  obtained  by  using  Eq.  (1).  The  semi-numerical  model  combines  the  analytic  approach  of  the 
Saturn  model  with  minimal  input  from  the  boundary  element  model,  thereby  using  analytic  current 
vs  gate  voltage  formulae  from  the  Saturn  1(V)  model.  The  angle-dependent  field  along  the  surface 
provided  by  a  geometric  parameterization  table  constructed  from  the  boundary  element  method. 
The  resulting  simulation  tool  is  then  in  a  position  to  rapidly  form  estimations  of  and  as  well 
as  other  parameters  (such  as  capacitance)  needed  in  a  full-scale  simulation  of  total  inductive  output 
amplifier  (lOA)  device  performance.  The  full-scale  device  simulation  is  a  program  of  considerable 
activity  in  the  FEA-RF  Amplifier  program  at  NRL  because  of  its  utility  in  relating  device 
performance  to  tip  parameters  to  affect  the  selection  process  for  FEAs  fabricated  by  widely  different 
techniques  [4]. 

FIGURES  OF  MERIT  AND  THE  FEA  UNIT  CELL 

Before  delving  into  the  theory  and  simulation  of  the  unit  cell,  the  practical  objectives  should  be 
identified.  The  performance  required  of  an  FEA  cathode  derives  from  the  specifications  of  the 
amplifier  for  which  it  is  intended.  In  the  case  of  an  inductive  output  amplifier,  which  requires  a  high- 
frequency  gated  cathode,  knowledge  of  the  I(V)  characteristic  curve  and  the  equivalent  circuit  of  a 
unit  cell  is  sufficient  to  predict  the  FEA's  effect  on  the  performance  of  the  amplifier.  Moreover,  the 
critical  information  is  conveyed  by  just  four  parameters:  Af„  and  By„,  the  maximum  available  current 
per  tip,  ipk,  and  the  net  gate-to-base  capacitance  per  unit  cell  (C).  General  guidance  for  optimizing 
these  parameters  can  be  obtained  simply  by  minimizing  the  drive  power  required  to  gate  the  beam. 

The  transconductance  of  a  voltage-controlled  current  source  is  defined  as  the  derivative  of  beam 
current  with  gate  potential,  i.e.,  =  3/beam/3'^g-  To  relate  FEA  performance  to  amplifier  gain,  a 

power  transconductance  is  defined  as  the  incremental  current  for  an  increment  in  drive  power:  gp  = 
3/beam/^^4r  amperes  per  watt).  In  a  gated  FEA  cathode,  the  gate  and  the  emitting  surface 

form  a  largely  capacitive  load  on  the  input  circuit;  the  power  required  for  a  voltage  swing  the 
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gate  is  just  -  (V2Q)(oCV;.y2,  where  C  is  the  capacitance  of  the  array,  and  the  voltage  swing  is 
defined  by  Vg  =  +  V^sinCoor),  in  which  (B  is  the  drive  frequency.  Consequently, 


dib 


2Q 


dr 


(12a) 


where  Q  is  the  quality  factor  of  the  input  circuit  and  ©o  is  the  resonance  frequency  of  the  input 
circuit.  (A  perfect  and  lossless  impedance  match  is  presumed.)  We  proceed  by  calculating  g„  at  the 
maximum  of  the  emission;  this  is  the  most  reasonable  approximation  for  strongly  modulated  beams 
because  the  input  power  requirement  is  driven  by  the  need  to  reach  that  maximum  emission.  By 
assuming  that  /beam  is  given  by  Eq.  (1 1)  and  letting  the  peak  value  of  /beam  Vg  be  denoted  by  Ip^ 
and  Vph  respectively,  we  obtain 


g  _  ^ovg  ^  ^ftt  ^ _ 2  Q _ 

{javg^  ^  ^  pk  ^tips^ti^ 


(12b) 


where  Ndps  is  the  total  number  of  unit  cells,  and  Cadd  is  capacitance  of  the  FEA  over  and  above  the 
sum  over  the  unit  cells.  The  rf  voltage  can  be  obtained  in  terms  of  the  modulation  current  ratio 
lavg/lpk .  by  convoluting  of  the  gate  voltage  signal  on  Eq.  (11).  For  lavg^tpk  <  0-4,  the  result  is: 


yielding  finally: 


'^^[la.gHpk] 


pk 


8p 


4ftQAfi,exp[-Bfi,/VpA 


2-h 


B 


Jh 


V 


(12c) 


(13) 


This  relation  offers  several  hints  for  the  optimization  of  FEAs  for  maximum  gain.  The  I(V) 
characteristic  curve  and  the  capacitance  are  the  FEA  characteristics  that  appear  here;  in  other  words, 
the  FEA  is  effectively  described  by  A/„,  Bf„,  and  C,  with  the  operating  point  given  either  by  ipk  or 
Vpk-  The  power  transconductance  is  most  sensitive  to  the  exponential  implicit  in  Eq.  (11),  and 
reducing  Bf„  is  the  most  effective  means  of  improving  the  FEA.  The  two  parameters  A/„  and  C 
occur  only  as  the  ratio  Afn/  C;  they  can  be  traded-off  against  each  other.  This  useful  figure  of  merit 
for  the  emission  gating  assembly  correlates  closely  to  the  gain  of  the  amplifier. 

ANALYTIC  EVALUATION  OF  CURRENT  VS  GATE  VOLTAGE  AND  FEA  PARAMETERS 

The  analytical  dependence  of  Af„  and  Bf„  on  material  and  geometric  parameters  for  simple 
geometries  must  be  determined.  (The  technique  is  extensible  to  more  complicated  surfaces,  but  this 
is  deferred  to  a  future  work.)  The  dependencies  rely  on  the  model  of  the  triode  geometry  (elliptical, 
hyperbolic,  sphere  on  post);  in  this  section,  we  concentrate  on  features  common  to  all  models. 
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Extraction  of  /(V)  from  J{F) 

The  total  current  from  an  arbitrarily  shaped  emitter  surface  is  given  by 

(14) 


where  J  is  the  current  density  and  dS  is  the  differential  surface  area.  A  widely  used  method  [20-23] 
of  estimating  total  current  from  a  tip  is  to  approximate  the  component  of  the  current  density  J 
normal  to  the  surface  by  the  1 -dimensional  Fowler-Nordheim  J{F)  as  given  in  Eq.  (9)  with  the 
approximations  for  v(y)  and  t(y)  as  given  in  Table  1.  We  are  more  careful  than  that  here,  in  that  we 
calculate  v(y)  and  r(y)  exactly. 

If  we  further  assume  rotational  symmetry  (the  cylindrical  coordinates  will  therefore  be  designated 
by  (p.z))  about  the  z  axis,  then  we  have 


ra. 


KVg)  =  2n 


1  + 


j[F{9)\dp, 


(15a) 


where  Zg(p)  describes  the  surface  of  the  emitter.  For  an  elliptical  emitter  (which  serves  as  a  generic  tip 
configuration  in  that  it  can  be  matched  to  different  models)  characterized  by  p  =  sin(9)  and  z  =  R 
Ug  cos(0),  where  is  the  base  radius  of  the  emitter,  this  becomes 


I{Vg)  =  271  al  JJ  +  y[F(x)]  dx , 


(15b) 


where  x  =  cos(0)  [27].  Hyperbolic  and  sphere  on  cone  models  result  in  analogous  formulae. 
Extension  of  the  FN  Equation  to  Semiconductor  or  Metallic  Multidimensional  Tips 

Three  complications  arise  in  the  treatment  of  real  emitters: 

•  Equation  (9)  as  given  ignores  temperature  and  variations  in  chemical  potential,  which  can  lead 
to  errors  when  applied  to  semiconductors  (although  not  for  metals). 

•  For  semiconductors,  or  for  metals  treated  with,  e.g.,  cesium,  the  barrier  height  experienced  by 
the  tunneling  electrons  may  be  so  low  at  high  fields  that  the  Fowler-Nordheim  equation  is  no 
longer  appropriate. 

•  Curvature  of  the  emitter  surface  alters  the  form  of  the  image  charge  from  Qlx  to  something 
that  is  geometry-dependent. 

•  Actual  tips  have  bumps  or  protrusions;  in  the  case  of  semiconductor  tips,  the  geometry  can  be 
pyramidal  rather  than  conical,  generating  a  field-enhancement  effect  along  the  pyramid  edge. 

For  the  first  two,  modified  Fowler-Nordheim  equations  can  be  found  for  both  metallic  and 
semiconductor  parameters  [19],  but  for  low  barrier  heights,  a  different  approach  toward  quantum 
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Table  1  —  Constants  and  Parameters 


Symbol 

Meaning 

Value  /  Equation 

Units 

Units 

A 

angstroms 

10-10 

meters 

eV 

electron  volts 

1.6x10-19 

Joules 

fi 

femtoseconds 

10-15 

seconds 

e 

unit  charge 

1.6x10-19 

Coulombs 

elifs  A^) 

current 

1.6x1012 

amps/cm2 

Fundamental  Constants 

c 

velocity  of  light 

2998.79 

A/fs 

a 

fine  structure  constant 

1/137.04 

- 

n 

Plank's  constant  /  2n 

0.658028 

eVfs 

mo 

electron  rest  mass 

511000 

eV/c2 

ao 

Bohr  radius 

0.529 

A 

Ry 

Rydberg  energy 

13.606 

eV 

^0 

permittivity  of  free  space  =  e^/47tafic 

1/180.95 

eV  A 

Statistical  Models 

A,- 

FN  parameter  of  ith  emitter 

- 

MA/V2 

Bi 

FN  parameter  of  ith  emitter 

- 

volts 

</> 

average  current  per  tip 

- 

MA 

N 

number  of  emitters 

- 

- 

Nb 

number  of  different  B  values 

- 

- 

ni 

number  of  emitters  with  B  =  Bi 

- 

- 

<T 

Gaussian  distribution  parameter 

- 

volts 

Bo 

reference  B  value 

- 

volts 

x± 

scaled  maximum  /  minimum  B  value 

(B±^Bo)/<T 

- 

Xo 

scaled  reference  B  value 

(x+  +  x.)/2 

- 

A 

scaled  spread 

(x+  -  x.)/2 

- 

V 

scaled  Gaussian  distrib.  parameter 

a/(2V) 

- 

Pk(z) 

Legendre  polynomial 

- 

Ck 

Legendre  polynomial  coefficient 

- 

- 
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Table  1  (Continues)  —  Constants  and  Parameters 


Fowler-Nordheim  Terms 


ypN 

Fowler-Nordheim  current  density 

OA,  exp(-bfi/F) 

Amps/cm2 

^fn 

parameter  (F  dependent) 

Nc 

Amps/cm2 

16  n^fi(pt{y)^ 

bfn 

J  ffsj{F)  parameter  (F  dependent) 

4  1 
3  1 

\  J 

eV/A 

v(y) 

function,  approximated  by  quadratic 

- 

^0 

constant 

=  0.93685 

- 

t(y) 

function,  approximated  by  constant 

=  1.05657 

- 

F 

field 

- 

eV/A 

<(> 

work  function  (metal) 

- 

eV 

chemical  potential  or  Fermi  level  (metal) 

- 

eV 

X 

electron  affinity  (semiconductor) 

- 

eV 

y 

Fowler-Nordheim  variable 

V(4QF)/«}) 

- 

Q 

planar  image  charge  coefficient 

ahc  1 

eV  A 

4  K,+  \ 

Ks 

dielectric  constant 

oo  for  metals 

- 

Nc 

conversion  factor 

1.6  X  10l2 

(amps/cm^)  x 

(e/A^  fs)’^ 

FEA  Device  and  Geometrical  Factors 

/(V 

current  from  unit  cell 

^fn 

V/exp[-FyV^] 

amps 

^fn 

FN  current  parameter  (constant) 

amps/cm^ 

^fn 

current  parameter  (constant) 

- 

volts 

gate  voltage 

- 

volts 

anode  voltage 

- 

volts 

^tip 

field  at  apex  of  tip 

F(a„0) 

eV/A 

X 

field  variation  factor 

- 

« 

a(0) 

radius  of  curvature  along  emitter 

- 

A 

b 

area  factor 

/(vp/y(Ftip) 

cm^ 

P 

inverse  temperature 

11604.5/r 

1/eV 

T 

temperature 

- 

Kelvin 

“V 

field  enhancement  factor 

- 

1/A 

c 

capacitance 

- 

Coul/volts 

h 

peak  current 

- 

amps 

^ave 

average  current 

- 

amps 

# 
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Table  1  (Continues)  —  Constants  and  Parameters 


Saturn  Model 


Fa 

field  coefficient  due  to  anode 

ValZa 

field  coefficient  due  to  gate 

Q^r^s 

^8 

charge  on  equivalent  gate  ring 

- 

r 

spherical  coordinate 

- 

e 

spherical  coordinate 

- 

p 

cylindrical  coordinate 

rcos0 

z 

cylindrical  coordinate 

r  sin9 

radius  of  gate 

- 

radius  of  ring  of  charge 

- 

^g 

radial  sphere  to  ring  separation 

+  a^2) 

radial  sphere  to  gate  separation 

^(Zg^  +  a/) 

cos(a) 

ring  angle 

t 

“thickness”  of  gate 

Og 

radius  of  equivalent  sphere 

- 

h 

distance  from  equivalent  sphere  to  gate 
plane 

' 

Za 

distance  from  equivalent  sphere  to  anode 

- 

Kip) 

complete  elliptical  integral  of  1st  kind 

rn/l 

dx 

0  s/ 1  -  phin^  [x] 

P 

K(p)  parameter 

(4pp')^^  /y 

Y 

ring  potential  parameter 

[(p-hp')^  +  (z-z')^]^^ 

eV/A 
eV/A 
eV  A 
A 

radians 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


A 
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transport  may  have  to  be  used  (Appendix  B).  With  regard  to  the  image  charge,  the  field  F  and  the 
barrier  height  %  can  be  altered  by  small  quantities,  dependent  on  the  radius  of  curvature  of  the  emitter 
at  the  location  where  the  current  density  J  is  being  evaluated  [27,  28].  To  a  very  good 
approximation,  %  and  F  can  be  replaced  by  the  quantities  F^  and  Xa  in  the  (modified)  Fowler- 
Nordheim  equation,  where 


F^(9)  =  F(e)  + 


[2a(e)]^ 


(16) 


so  that  the  barrier  potential  resembles  x^(e)  -  F^{Q)  x  -  Qlx.  In  the  models  considered  here,  the 
emitter  tip  is  approximated  by  a  sphere,  so  that  0(9)  =  flsphere  for  the  values  of  9  where  7(^(9))  is 
significant.  Finally,  with  regard  to  surface  curvature,  the  bumps  and  protrusions  on  axis  have  the 
dominant  effect.  Ridges,  being  a  protrusion  in  a  two-dimensional  space,  have  less  field  enhancement 
than  bumps  in  a  three-dimensional  space,  and  therefore,  pyramidal  complications  are  not 
immediately  of  concern  (Appendix  C). 

Fowler-Nordheim  Parametrization  of  I{Vg)  and  Area  Factor 

In  pursuit  of  an  analytic  I{Vg),  two  assumptions  are  expeditious:  Assume  that  the  current  density 
can  be  approximated  by  Eq.  (9),  (note  that  this  does  not  imply  a  simple  correspondence  with  Eq. 
(11),  as  v(y)  and  f(y)  are,  in  general,  field-dependent).  Also  assume,  at  least  for  small  9,  the  field  to  be 
approximated  by 


F(9)  = 


1  +  X  ( 1  -  cos  (9)]  ’ 


(17) 


where  —  field  on-axis  and  X  characterizes  the  falloff  of  field  along  the  equivalent  sphere 
representing  the  surface,  both  of  which  are  gate-voltage  dependent.  Inserting  these  approximations 
into  Eq.  (15b)  gives 


n 


■  2jta^  Jpjy[(F^p 


VF^-(/?^-l)(l-y)^ 


MO) 

Ptip 


Xy\dy 


(18) 


where  y  =  7-cos( 6)  and  a^  is  the  base  radius  of  the  elliptical  emitter.  The  coefficient  of  7fn(^o-  )  ‘s 
the  area  factor  (b).  Finally,  F„^,  X,  C,  6,  and  the  peak  and  average  currents  Ip  and  are  dependent 
on  the  geometrical  factors,  and  are,  therefore,  dependent  on  the  model  of  the  triode  we  adopt.  This  is 
taken  up  next. 

THE  SATURN  MODEL 

Potential  Distribution  of  Sphere  and  Ring,  and  the  Electric  Field 

In  some  sense,  the  simplest  model  of  an  ungated  field  emitter  (diode)  is  given  by  the  “floating 
sphere”  model  [24,  29],  in  which  a  sphere  whose  radius  is  equal  to  the  tip  radius  of  the  field  emitter 
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is  suspended  between  parallel  plates.  The  simplest  modification  to  this  model,  which  endeavors  to 
approximate  a  triode,  is  to  insert  a  ring  of  charge  (not  necessarily  coplanar)  somewhere  along  the 
symmetry  axis.  We  further  replace  the  parallel  plates  by  a  background  field  of  the  form  -Fz.  Let  the 
origin  of  the  coordinate  system  be  centered  on  the  sphere  approximating  the  tip  (Fig.  1).  The 
potential  anywhere  is  therefore  the  sum  of  three  components  —  the  anode,  the  gate  ring  (or  simply 
gate),  and  the  sphere. 


r  cos(  6) 


Fig,  1  — Saturn  model,  showing  relevant  parameters 


For  the  anode,  the  potential  is  given  by 

^a  =  -^a’'  COS  (0)  . 


(19) 


where  we  let  =  VJZa,  where  is  the  potential  of  the  anode  and  is  the  distance  from  the  center  of 
the  sphere  to  the  anode  plane.  The  justification  for  this  assignment  goes  back  to  the  diode  calcula¬ 
tions  in  which  an  analogous  expression  was  found  on  axis  [27].  Analogously,  and  using  the 
nomenclature  introduced  in  Fig.  1  and  Table  1,  we  define  the  gate  potential  to  be 


Qg  , 

sJo 

p2  +{z-Zg)^  +  a^  -  2pa^cos((p) 

2n  'fW 

71  Y 

O 

^  z 

rg  1  =  0 

1  P/(cos  (a))  P,(cos  (6))  , 

(20a) 


where  we  have  introduced  the  integral  definition,  the  complete  elliptical  integral  K{p),  and  the 
Legendre  polynomial  expansion  for  r  <  r^.  The  value  for  is,  as  yet,  unspecified;  note  that  it 
implicitly  includes  a  factor  of  (IMnCo).  K(p)  is  defined  as 
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Kip)  = 


J-n/2 

0 


_ dx _ 

\l \  -  p^  sin^jc 


(20b) 


Finally,  for  the  tip,  the  potential  is  given  by 

(f)j  =  A ;  r-  P;(cos  (9)) , 


(21) 


where  the  A;’s  are  to  be  determined. 


Determination  of  the  Legendre  Coemdents,  Gate  Charge,  and  Field  on  the  Sphere 

On  the  surface  of  the  sphere,  the  potential  (|)  =  +  (j)  +  (j)^  vanishes.  This  serves  to  specify  the 

A['s: 


QilaM‘+^)  (a  V  r  , 

^l  =  -Tj{-r)  (?j)  ^/H(a))-F^-^8/i  , 


(22) 


where  5//.  is  the  Kronecker  delta  function.  Next,  around  the  gate  ring,  the  equipotential  lines  are 
approximately  toroidal.  Let  a  point  separated  from  the  gate  ring  by  a  distance  r  =  be  at  the 

potential  Vg.  Assuming  t  is  not  large  compared  to  other  parameters,  it  can  be  taken  as  half  of  the  gate 
thickness  (this  approximation  is  poor  for  typical  FEAs  since  the  gate  thickness  is  not  in  general  small 
by  comparison  to  other  terms).  Then  Qg  is  determined  by 


y  -2.  Q  +F  z 


8 


Jo 


a 


1-1^ 

g 


‘  =  A'‘'gJ 


fV 


(23a) 


Note  that  r'g  =  V(Zg2  +  ^^2)  ^  +  ^2)  _  «  Og,  only  the  first  term  of  the  summation 

need  be  retained,  and  we  have  the  approximation 


iV  -F  z' 

2K(p,)  a, 


(23b) 


where  and  are  p  and  y  evaluated  at  p  =  and  z  =  Zg,  respectively.  The  charge  on  the  sphere  is 
given  by  A;  for  /  =  0,  or  =  Qg  /  r^.  As  expected,  this  result  is  analogous  to  the  image  charge 
solution  to  the  problem  of  a  conducting  sphere  under  the  influence  a  point  charge  [30]. 

Along  the  surface  of  the  sphere,  the  force  {F  =  -eE)  on  the  electrons  points  radially  outward.  We 
can  therefore  calculate  the  field  along  the  sphere  by  taking  the  radial  derivative  of  ([)  and  evaluating  it 
at  r  =  a^.  Doing  so  yields 
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F(a,.e)  S  -F(e) 

=  -3F^cos(e)-F  £  (2Z+1) 
*  /  =  0 


—  j  Pi  [cos  (a))  Pi  (cos  (6))  , 

Vs) 


(24) 


where  we  have  introduced  ,Fg  =  Qgl(rga^).  As  it  should,  Eq.  (24)  naturally  separates  into  a  cathode 
and  anode  part.  In  general,  F„p  6  but  is  something  more  complicated.  For  a  distant  anode 
generating  a  relatively  weak  F„,  the  approximation  acquires  merit.  The  field  at  the  tip  (0  =  0°)  can  be 
found  by  replacing  F;(cos(0))  with  (1)  in  Eq.  (24).  Again,  if  we  assume  that  the  sphere  radius  is 
much  smaller  than  the  gate  radius,  then  we  can  approximate  the  X  factor  by 


X,  =  3 


F  +F 


^  "*■  ^g 


a,Zg' 


JPupzfl 

^tip 


(25) 


where  the  s  subscript  on  X  refers  to  Saturn.  One  consequence  of  Eq.  (24)  is  that  the  field  is  not 
maximized  for  the  tip  of  the  emitter  lining  up  with  the  gate  plane.  The  emitter  should  be  below  (or 
above)  the  gate  plane  (Zg  >  a^)  by  an  amount  that  maximizes  F  as  a  function  of  Zg-  This  has  been 
previously  noted  in  numerical  experiments  [20].  Finally,  note  that  A5  will  be  small  if  the  quantity 
{asZg/fg^)  is  small,  as  will  happen  in  practice  when  real  FEAs  are  considered  for  which  the  numerator 
is  vanishing  by  comparison  to  the  denominator. 

Current  and  Capacitance  in  the  Saturn  Model 


For  the  Saturn  model,  /?  =  1,  for  which  Eq.  (18)  reduces  to  (dropping  the  s  on  X): 


ri 


I(Vg)  =  2na^Jj^Fiip) 


0  (l+M 


2  exp 


bjniO) 


Xy 


dy. 


(26) 


where  bf„(0}  is  bf„  evaluated  zt  F  =  0.  A  good  approximation  to  I(Vg)  is  then  given  by 


/(V  =  ^F2[x,%^]7fn(F„P  ’ 

A.  ^  t^tip 


where  F„(k^)  is  given  by  [27]; 


F„{Kx)=j^(l+y)-"e-^dy  . 


(27a) 


(27b) 


F2(X,x)  is  bounded  from  below  by  l/x  (small  field  limit)  and  from  above  by  X/(l+X)  (large  field 
limit).  For  small  fields  at  the  tip,  we  therefore  see  that,  in  contradiction  to  the  current  density  J,  for 
the  current  /,  \n{I/F^)  is  approximately  linear  in  1/F.  The  coefficient  of  Jfn(F)  in  Eq.  (27a)  is 
typically  referred  to  as  the  “area  factor.” 

The  capacitance  in  the  Saturn  model  is  given  by  the  ratio  of  the  charge  on  the  sphere  and  the 
gate  voltage  Vg,  which  is  related  to  the  charge  on  the  ring  by  the  term  in  the  Legendre  expansion. 
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The  charge  on  the  ring  is  given  by  Eq.  (23b).  An  analytic  expression  exists  for  the  sphere-to-gate 
capacitance  that  can  be  found  in  a  manner  analogous  to  that  used  for  the  boundary  element  model 
(described  below).  We  note  here  that,  depending  on  the  geometry,  Eq.  (27)  and  the  expression  for  C, 
may  or  may  not  be  adequately  approximated  by  the  Saturn  model. 

Two  points  worth  mentioning  regarding  Eq.  (27).  First,  Eq.  (27a)  owes  its  simple  form  to  the 
approximation  v(y)  -  used  in  Eq.  (26).  Two  problems  are  that  v(y)  is  not  linear  in  Fup,  nor  is 

Vo  -  y2  the  best  (least-squares)  approximation  to  v(y)  (Appendix  B).  The  former  problem  is  resolved 
by  noting  that  for  field  values  of  interest,  a  linear  approximation  to  v(y)  in  the  field  is  a  reasonable 
approximation.  For  the  latter,  v(y)  may  always  be  Taylor  expanded  about  a  (typical)  field  of  interest 
in  order  to  generate  v(y)  =  v  ^  -  v;  y^,  thereby  entailing  a  minor  modification  to  Eq.  (27a).  Second, 
JpNiPtip)  in  Eq.  (27a)  can  be  replaced  by  any  other  (superior)  quantum  transport  model,  e.g.,  one 
based  on  the  Airy  function  solution  to  Schrodinger's  equation  [33].  This  is  particularly  important  if 
the  work  function  becomes  so  low  that  at  high  fields,  the  Fowler-Nordheim  equation  is  no  longer  a 
good  approximation  to  J(F).  We  consider  a  case  of  this  below  when  the  experimental  data  from 
MIT-Lincoln  Labs  are  discussed. 

Deficiencies  of  the  Saturn  Model 

For  a  given  F tip  and  X,  we  have  found  that  the  Saturn  model,  i.e.,  Eq.  (27),  is  a  good 
approximation  to  Eq.  (26).  The  problem  with  the  Saturn  model  lies  with  its  estimate  of  Fjip,  which 
ultimately  depends  on  the  estimate  of  Q^,  and  also  its  estimation  of  X.  Although  it  is  true  that  the 
charge  density  spikes  near  an  edge  [30]  giving  the  justification  for  replacing  the  gate  (a  sheet  of 
charge  with  a  hole  excised  from  it)  with  a  ring  of  charge  in  the  proximity  of  the  edge  of  the  hole,  the 
remainder  of  the  gate  (and  the  charges  on  it)  cannot  be  neglected.  Consequently,  for  a  given  the 
Saturn  model  tends  to  overestimate  the  value  of  Vg  associated  with  it,  as  the  result  of  to  its  neglect  of 
the  remainder  of  the  gate  compared  to  a  full  numerical  solution.  Attempting  to  correct  this  limitation 
results  in  the  addition  of  more  rings  of  greater  radii  with  different  Q’s,  but  this  is  too  similar  to  the 
actual  numerical  solution  to  merit  elaboration.  In  summary,  once  a  (correct)  estimate  between 
and  Vg  is  established,  the  remainder  of  the  machinery  adopted  for  the  Saturn  model  can  be  used  to 
estimate  current  and  related  quantities. 

THE  BOUNDARY  ELEMENT  MODEL 

Geometry,  Parameters,  and  Implementation 

A  model  that  corresponds  well  to  the  actual  geometry  of  the  gated  field  emitter  is  expected  to 
provide  a  more  accurate  comparison  with  the  experimentally  measured  quantities.  The  axially 
symmetric  unit  cell  model  selected  consists  of  an  anode,  a  gate  with  hole,  and  a  base  plane  with  an 
emitter  tip  protruding  (Fig.  2).  This  model  is  representative  of  several  of  the  "vertical  emitter" 
structures  currently  being  considered  in  the  FEA-based  RF  Amplifier  Program  at  NRL.  The  small 
sizes  characteristic  of  the  apex  of  the  emitter  tips  in  comparison  to  the  size  and  distances  of  the  other 
electrodes  suggest  the  use  of  a  nonuniform  discretization  of  the  computational  domain.  The  use  of 
typical  finite-difference  formulations  with  uniform  mesh  is  therefore  contraindicated  (as  in  SIMION 
and  EGUN  simulations).  A  variable-mesh  finite-difference  or  finite-element  technique  could  be  used 
[21-23].  However,  because  the  principal  item  of  interest  is  the  emission  at  the  boundary,  particularly 
in  the  small  region  near  the  emitter  tip,  a  boundary-element  model  was  chosen. 
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Fig.  2  —  Generic  FEA  unit  cell,  showing  relevant  parameters 


The  boundary  element  technique  [24]  discretizes  the  boundary  into  sections  (in  our  model, 
annular  ribbons).  These  sections  are  taken  to  be  quite  small  in  the  vicinity  of  edges  and  comers  and 
increase  in  size  (quadratically,  in  the  example  to  be  shown  later)  in  the  smooth  boundary  areas.  The 
list  of  boundary  elements  and  their  physical  attributes  is  developed  from  a  parameterization 
suggested  by  experimentally  pertinent  quantities.  These  qualities  include  work  function,  potential,  tip 
radius,  tip  height,  gate  hole  radius,  base-gate  distance,  gate-anode  distance,  gate  thickness,  and  type  of 
tip  (sphere  on  cone,  ellipsoid,  tip  on  post,  post  height).  A  graphical  user  interface  allows  a  rapid 
analysis  of  the  system  because  of  integration  of  the  boundary  element  generator,  potential  solver, 
semi-numerical  post-processor,  and  plotting  package.  After  the  geometry  is  specified,  the  time  taken 
by  the  program  package  to  extract  the  Fowler-Nordheim  and  as  well  as  Cpgyj^,  parameters  is 
less  than  a  minute  on  a  desktop  computer. 

Calculation  of  the  Potential 


The  potential  anywhere  can  be  calculated  by  integrating  over  the  surface  charge  densities  on 
electrodes: 


(28a) 


where  dO.  is  the  differential  surface  element.  Consider  the  approximation  in  which  the  conducting 
boundaries  are  broken  up  into  ribbons  of  charge,  each  with  a  constant  surface  charge  density  [25]. 
This  tends  to  lead  to  a  faster  numerical  solution  than  assuming  that  the  surface  charge  density  has  a 
linear  dependence  [26]  and  reflects  our  eventual  intent  to  rely  on  a  hybrid  between  the  Saturn  model 
and  the  boundary  element  model.  Consequently,  Eq.  (28a)  assumes  the  discretized  form: 
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where  ^  -  P/),  cr/^1/2  is  the  charge  density  on  the  iih  ribbon,  p  =  V(4pp')/Y  and  y 

=  V[(p+p')2  +  (z-zO^].  The  p,-  are  quadratically  spaced,  so  that  the  ribbons  are  thinner  near  the  axis 
of  symmetry  or  near  a  conductor  edge. 


At  the  midpoint  of  each  ribbon,  (f)  is  specified  by  the  potential  of  the  conducting  surface.  There 
are  N  ribbons,  and  we  therefore  have  N  values  of  4»  denoted  by  [(])],•  =  (t)(p,+i/2,  z,+i/2).  where  (j)  is  an  N- 
component  vector,  p,+i/2  =  p,-  +  e,/2,  z,+i/2  =  z,-  +  5,+i/2e,/2,  and  e,-  =  (p,+i  -  p,).  We  can  therefore 
rewrite  Eq.  (28b)  as  (where,  analogous  to  the  potential  vector,  [a],-  =  di+m) 


[M]  ij  -  1  12 


p  Y 


(29) 


where  p,-  and  z,-  are  implicit  in  the  definition  of  p  and  y.  The  (ji-vector  is  a  known  quantity  (it  is  the 
value  of  the  voltage  segments  on  the  tip,  gate,  and  anode),  and  can  be  calculated.  Eq.  (29)  can 
then  be  inverted  to  give  the  values  o^the  charge  densities  on  each  ribbon:  a  =  M-'  •  ()).  An 
advantage  of  this  formulation  is  that  M  need  be  calculated  only  once  (being  geometry,  and  not 
potential,  dependent)  and  inverted  once  (using,  for  example,  LU  decomposition);  thereafter,  for  a 
different  voltage  on  the  gate  and  anode,  only  the  (j)  vector  needs  to  be  altered  to  obtain  a  new  ct. 


The  integral  in  Eq.  (29)  is  generally  evaluated  by  using  a  4-point  Gaussian  quadrature  routine, 
except  when  li-yl  <  3,  in  which  case  a  10-point  Gaussian  quadrature  routine  is  used.  The  ring 
potential  of  Eq.  (20a)  can  be  derived  from  Eq.  (28)  by  evaluating  the  integrand  at  the  midpoint. 
When  /  =7,  K(p)  contains  a  logarithmic  singularity  of  the  form  (1/2)  ln[16/(l-p2)]  as  p  approaches  1. 
The  singular  portion  of  the  integrand  can  be  analytically  integrated.  Dropping  the  (/+1/2)  subscript 
on  (j)  and  (p)  and  the  (/)  subscript  on  (e)  on  the  right-hand  side,  we  have: 


P-iln 

7  2*" 
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1  -p‘ 


dp'  = 


288p^ 


^[96p2-(j2  +  2)e2] 


3  In 


^  256  p2  ^ 

^(^2+  i)e2^ 


+  2 


+  96p2 

i 


(30) 


The  remainder,  i.e.,  p'  {Kip)  —  (1/2)  ln[16/(l— p2)]}/Y^  can  be  integrated  by  Gaussian  quadrature. 

Once  the  surface  charge  densities  a, +1/2  are  found  by  inverting  the  matrix  equation  (Eq.  (29)), 
the  field  at  the  surface  F,-  and  (normalized)  charge  per  ribbon  2,  can  be  evaluated: 


and 


~  2 *  ^2  ^  1/2  P'+  1  /2 (32) 

(note  the  distinction  between  the  grid  spacing  e,-  and  the  permittivity  of  free  space  e^).  The  field  F,-  is 
constant  along  the  ith  ribbon.  From  the  F,’s,  the  total  current  from  the  emitter  can  be  found  by 
summing  over  the  current  contributions  from  each  ribbon  that  constitute  the  tip: 
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tip 


/(V  )  =  2 71  y  P/+ 1/2 £/ m  ■ 


(33) 


Finally,  if  the  potential  away  from  the  electrode  surfaces  is  needed,  as  for  field  lines  and  trajectories 
[13],  then  a  useful  approximation  is  to  replace  each  ribbon  by  a  ring  of  charge  located  at  p,+i/2.  with 
a  total  charge  g,-.  The  potential  at  any  point  (p,z)  can  then  be  calculated  by: 


N 

(t>(p,z)  =  X, 

I  =  1 


27t 


r2n 


</cp 


s/p^  +  (z- Zi+ I, if +  pj+ 1/2  -  2pP,  +  l/2COS((p) 


(34) 


o  2^(A) 
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where,  as  before,  the  e^-dependent  factors  are  implicit  in  the  Q’s.  Finally,  the  capacitances  are  found 
from  the  matrix  equation  [(t,b,g,a)  refers  to  (tip.  base,  gate,  anode),  respectively]: 


Qt 

^tt  ^tb  ^tg  ^ta 

‘t'r 

Qb 

^bt  ^bb  ^bb  ^bb 

• 

Qg 

^gt  ^gb  ^gg  ^ga 

Qa 

^at  ^ab  ^ag  ^aa 

<!>« 

where  Cjj  =  Cj^,  although  this  is  not  shown  explicitly.  To  evaluate  the  C’s,  all  the  potentials  ((]),•) 
except  one  (^j)  are  variously  set  equal  to  0;  we  then  have  C,y  =  The  charges  are  found  by 

summing  over  the  ribbon  charges  (as  in  Eq.  (32)),  e.g.,  for  the  tip,  gfip  =  Qi  fo/"  ^  ^*P}- 
alternate  approach,  which  is  used  here,  is  to  calculate  C(Vg)  =  C,g  +  CfVJVg,  Cfg  will  then  be  the  y- 
intercept  in  a  plot  of  C  vs  l/V^,  and  the  tip-to-gate  capacitance  is  obtained  without  separately 
calculating  the  capacitance  matrix. 

THE  SEMI-NUMERICAL  MODEL 

In  the  semi-numerical  model,  we  desire  to  retain  the  simple  elegance  of  Eq.  (27)  for  the  /(Fg) 
relation.  From  the  qualitatively  correct  behavior  of  X.,  (Eq.  (17))  and  F,ip  (derivable  from  Eq.  (24)), 
we  are  in  a  position  to  estimate  the  parameters  that  will  most  dramatically  affect  the  I(V)  relations  (the 
behavior  becomes  more  quantitatively  correct  if  the  gate  thickness  parameter  t  is  very  small;  when  the 
Saturn  model’s  charge  ring  location  no  longer  matches  the  "mean"  radius  of  the  charge  distribution 
on  the  gate,  or  when  the  anode  influences  Qg,  the  approximation  fails).  In  this  section,  we  provide  a 
numerical  approach  for  obtaining  a  refined  estimate  of  the  X  and  parameters  by  using  results  of 
a  minimal  set  of  boundary  element  computations,  and  thereby  obtain  the  most  expeditious  means  of 
calculating  the  Af„  and  Bf„  parameters  and  the  tip-to-gate  capacitance.  In  particular,  we  require  an 
algorithm  to  calculate  F,ip(Vg)  and  XfVg)  from  the  boundary  element  model  for  use  in  Eq.  (27),  in 
place  of  the  values  calculated  by  the  Saturn  model.  From  the  linearity  of  the  field  vs  gate  voltage 
relations,  it  is  clear  that  for  the  field,  we  have 

^tip(V'g)  =  B^Vg  +  y> 


(36) 
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(in  keeping  with  the  field  enhancement  factor  convention,  6  is  the  slope).  Although  not  so  obvious 
analytically  (although  apparent  graphically),  we  can  also  use 

>^(v,)  =  (6;,/vp  +  yx.  (37) 

where  the  subscripts  on  6  and  Y  refer  to  the  X(Vg)  and  F(Vg)  slopes  and  intercepts,  respectively.  A 
similar  trick  can  be  used  to  find  Q/V^)  to  evaluate  C,g  using  the  Saturn  model,  but  this  is  deferred  to 
a  future  work. 


Because  of  the  linearity  of  Eqs.  (36)  and  (37)  in  Vg  and  its  inverse,  only  two  boundary  element 
calculations  are  required  to  determine  the  6  and  Y  parameters  for  a  given  geometry,  likewise  with  the 
calculation  of  Cfg.  A  table  of  these  values  for  various  geometries  can  then  be  constructed  from  which 
any  particular  geometry  can  be  extrapolated.  Ideally,  the  Fowler-Nordheim  A  and  B  parameters 
should  be  evaluated  by  performing  a  linear  fit  of  \n[I/Vg'^]  vs  MVg  in  the  regime  of  interest. 
However,  the  semi-numerical  model,  in  conjunction  with  the  Saturn  model,  can  be  used  to  infer 
parametric  dependencies  of  A  and  B  on  the  geometry  and  work  function.  For  example,  assume  that 
v(y)  =  Vo  -y2,  and  anode  effects  are  negligible,  i.e.,  0.  From  Eq.  (27),  we  therefore  identify 


^  !^a3/2 


(38a) 


which  hearkens  to  the  field  enhancement  factor  typically  used.  An  analytical  approximation  for  Aaj 

is  more  problematic,  due  to  the  dependence  of  F„(A,x)  on  the  gate  voltage  (that  Af^  must  be  Vg- 

dependent  has  been  noted  before  [16,20]).  For  typical  parameters,  we  have  found  that  F„  satisfies 
/  ^ 

(x+n)"^  [7-exp(-A(jt-i-/i))]  <  F„(A,jc)  <  A/(i+A),  the  boundaries  representing  the  lower  and  upper 

ranges  of  the  gate  voltage.  Thus,  order  of  magnitude  estimates  of  A are  bounded  from  below  by: 


M  2 


m<l> 


5  g 


V.exp 


16  ^ 


2  m 


(38b) 


for  the  small  Vg  case,  and  where  the  (/)  subscript  on  p  has  been  suppressed.  For  the  large  Vg  case: 


16^ 

-2 - exp  —Q  /  — 

r  7t(A-t-l)(l)h  Y  4,h2j 


(38c) 


when  the  gate  voltage  approaches  its  upper  limit.  As  indicated  in  Eq.  (38b),  a  plot  of  Ay„  vs  Fg  is 
approximately  linear  in  the  low  gate  voltage  regime.  Equation  (38)  is  given  for  pedagogical 
purposes  and  for  order-of-magnitude  estimates.  As  a  means  of  calculating  the  Fowler-Nordheim  fit 
parameters  in  Eq.  (11),  the  voltage  dependence  of  Af„  in  Eq.  (38)  alters  the  best  estimates  of  the  FN 
A  and  B  parameters  as  would  be  obtained  from  a  linear  fit  performed  on  Eq.  (27)  (Appendix  A)  in  a 
given  Vg  regime;  ignoring  anode  effects  can  cause  differences. 

Although  the  Saturn  model  can  give  poor  quantitative  estimates  of  the  field  enhancement  factor 
p,  the  qualitative  dependencies  can  be  inferred.  In  the  Saturn  model,  the  field  at  the  tip  primarily 
varies  to  first  order  as  the  ratio  of  the  charge  on  the  gate  ring  with  the  product  of  the  distance  from 
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the  tip  to  gate  ring  and  the  tip  radius,  i.e.,  F^p  ^Qg/a^rg^  pVg,  where  we  have  assumed  that  the  tip 
radius  is  much  smaller  than  the  ring  radius.  If  we  assume  that  the  tip  of  the  emitter  sphere  lies 
approximately  within  the  ring  plane,  that  the  anode  is  negligible,  and  that  the  asymptotic  form  of  Eq. 
(23b)  is  used,  we  find: 


^tip 


71 


a,  In 


(38d) 


At  this  point,  t  is  an  unknown  parameter  that  can  even  vary  with  gate  voltage,  although  it  is  related  in 
some  sense  to  the  gate  parameters.  The  coefficient  of  Vg  is  the  "Beta  factor"  ()8).  The  relationship 
between  the  field  at  the  tip  and  the  charge  on  the  gate  ring  allows  us  to  conclude  from  Eq.  (38d)  that 


na, 


(38e) 


To  reiterate,  Eq.  (38)  shows  the  nature  of  the  dependences  of  Af„,  Bf„,  ^tipy  and  Ctg  on  material  and 
geometry-dependent  parameters.  Generally,  neglecting  the  anode  and  complications  due  to  finite 
gate  width,  tip  cone  angle,  and  so  on,  causes  the  quantitative  estimates  based  on  the  above  equations  to 
be  not  reliable  except  for  particular  cases.  However,  as  qualitative  descriptions,  Eqs.  (38a-e)  are  in 
fact  useful;  they  suggest  the  manner  in  which  capacitance,  field,  and  so  on  vary  with  geometrical 
parameters. 


NUMERICAL  RESULTS 


Objectives 

In  this  section,  we  undertake  six  tasks: 

•  use  statistical  analysis  with  boundary  element  computations  to  examine  experimental  data  to 
develop  and  determine  the  applicability  of  the  methodology; 

•  identify  where  improvements  can  be  made  and  under  what  circumstances  the  conditions 
requiring  improvement  arise; 

•  demonstrate  that  the  analytic  generalizations  of  the  statistical  and  Saturn  sections  are 
predictive; 

•  define  the  protocol  behind  the  application  of  the  semi-numerical  method; 

•  demonstrate  the  utility  of  the  boundary  element  method  for  analyses  beyond  current-voltage 
characterization,  i.e.,  trajectory  determinations  and  field  lines  within  the  unit  cell;  and 

•  delineate  the  incorporation  of  the  semi-numerical  method  into  more  extensive  simulation 
programs. 

Protocol  for  the  Comparison  of  Theory  with  Experiment 

To  judge  the  utility  of  the  boundary  element  model  in  correlating  with  experimental  data,  we 
compare  the  results  of  the  boundary  element  simulation  based  on  the  geometry  shown  in  Fig.  2  with 
experimental  data  from  two  sources:  the  first  has  been  provided  by  the  Massachusetts  Institute  of 
Technology  Lincoln  Laboratories  [12]  (Case  1),  and  the  second  is  characteristic  of  the  SRI  "Low 
Frequency  Cathode"  discussed  in  the  literature  [9]  (Case  2).  The  low  gate  voltage  regime  of  the  Case 
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1  experimental  data  was  used  in  the  least-squares  fitting,  (it  has  been  argued  that  the  high  voltage 
regime  is  affected  by  space  charge). 

Several  effects  complicate  the  comparison  of  theory  with  experiment: 

•  the  actual  work  function  and  effective  tip  radius  can  be,  and  are,  affected  by  the  history  of 
operation  of  the  FEA  and  the  processing  steps,  e.g.,  “plasma  cleaning,”  and  contaminants 
present; 

•  the  actual  FEA  surface  is  not  smooth,  as  implied  by  the  theoretical  procedure — there  is  always 
the  possibility  of  microprotrusions.  The  scale  of  the  microprotrusions,  or  surface  non¬ 
uniformities,  will  introduce  another  parameter  into  the  theory. 

•  the  experimental  current  per  tip  is  obtained  by  taking  the  ratio  of  the  total  FEA  current  and  the 
total  number  of  tips  in  the  array. 

•  the  shape  of  the  gate  in  the  FEA  is  more  complex  than  assumed  in  the  model  here  (although 
this  can  be  accounted  for  in  the  future  by  altering  the  geometrical  parameterization); 

•  the  position  and  potential  of  the  “virtual  anode”  will  not,  in  general,  be  at  the  same  location 
and  potential  over  a  range  of  gate  voltages,  as  we  have  assumed  here  (this  can  be  accounted  for 
in  a  more  careful  treatment);  and 

•  the  distribution  of  emitters  in  all  probability  is  more  complex  than  the  simple  distributions 
examined  here. 

The  first  two  concerns  are  dealt  with  by  appealing  to  “effective  work  functions”  and  “effective 
tip  radii.”  Insofar  as  the  work  function  of  the  tip  is  affected  by  contaminants,  we  treat  it  as  an 
adjustable  parameter  (within  bounds).  Microprotrusions  and  nonspherical  tips  are  dealt  with 
analogously:  the  analysis  here  returns  a  length  scale  that  is  characteristic  of  the  tips  in  question,  not 
the  actual  tip  radius.  As  there  is  some  evidence  that  surface  non-uniformities  and  conditions  can 
enhance  emission,  we  might  expect  that  surface  protrusions  and  conditions  for  FEAs  made  by 
different  processes  may  nevertheless  result  in  comparable  effective  tip  radius  estimates. 

Regarding  the  position  and  potential  of  the  anode:  In  simulations  of  FEA  performance,  the 
equipotential  lines  at  a  distance  of  a  few  gate  hole  diameters  from  the  gate  plane  are  seen  to  be 
approximately  linear.  We  are  therefore  Justified  in  replacing  one  of  these  equipotential  lines  with  a 
“virtual  anode”  held  at  the  same  potential.  This  approach  is  of  great  utility  in  the  present 
simulations  due  to  the  finite  size  of  the  unit  cell  and  the  desire  to  keep  the  computational 
requirements  to  a  minimum,  that  is,  we  desire  the  simulations  to  run  on  a  desktop  computer. 
Furthermore,  it  is  justified  insofar  as  the  effect  of  the  anode  is  overshadowed  by  the  gate  (although  it 
is  not  insignificant).  However,  the  position  and  potential  of  the  actual  anode  is  usually  not  reported, 
or  is  uncertain,  making  the  potential  and  location  of  the  virtual  anode  a  problem.  In  Case  1,  the 
actual  anode  was  reported  to  have  been  held  at  100  V  at  a  distance  of  1(X)  (tm  from  the  gate,  from 
which  we  approximate  that  the  virtual  anode  should  be  held  at  0.5  V  at  a  distance  of  0.5  |im  from  the 
gate  (the  virtual  anode  is  3.33  gate  diameters  away  from  the  gate  plane).  In  Case  2,  the  potential  of 
the  anode  was  reported,  but  the  anode  was  at  an  angle  with  respect  to  the  gate  plane;  we  therefore 
made  a  heuristic  estimation  as  to  the  potential  and  distance  away  of  the  virtual  anode. 

To  correlate  theory  with  experiment,  we  perform  the  following  adjustments  to  the  geometry  and 
material  parameters  in  order  to  obtain  correlation  with  experiment: 

•  The  value  of  Bf^  is  estimated  first  by  adjusting  the  tip  radius  until  approximate  agreement  with 
experiment  is  achieved.  This  “agreement”  should  actually  be  below  the  experimentally 
observed  value  if  we  intend  to  augment  the  simulation  with  a  statistical  analyis. 
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•  Depending  on  the  statistical  model  used,  we  next  endeavor  to  match  the  Afn  parameter  by 
adjusting  the  distribution  of  tips.  Typically,  the  wider  the  distribution,  the  more  Afn  is 
decreased,  but  also,  the  more  Bfn  is  increased  (thereby  explaining  why  we  desire  to  undershoot 
the  Bfn  estimate  above).  The  introduction  of  the  statistical  analysis  parameters  primarily 
reduces  the  Afn  value,  and  therefore  is  the  primary  means  of  accounting  for  the  fact  that  the 
experimental  current  per  tip  is  always  smaller  than  the  current  characteristic  of  the  unit  cell 
(sometimes  significantly). 

•  The  statistical  analysis  generally  changes  the  Bfn  estimate.  We  therefore  return  to  the  first  step, 
and,  depending  on  circumstances,  tweak  the  work  function  and  the  virtual  anode  estimates  in 
addition  to  the  tip  radius,  until  correlation  is  obtained. 

This  procedure  will  not  produce  a  unique  Afn  and  Bfn  estimate,  being  ad  hoc  in  its 
implementation.  Nevertheless,  we  have  found  that  for  the  experimental  cases  we  have  simulated,  the 
allowable  range  of  values  for  effective  tip  radius,  work  function,  and  virtual  anode  are  rather 
constrained.  For  example,  for  the  three  distributions  considered  below,  decreasing  ())  to  decrease  Bfn 
results  in  Afn  being  made  larger.  The  validity  of  the  modeling  effort  is  enhanced  by  this  observation, 
as  well  as  the  observations  that  the  work  function  converged  upon  is  within  the  range  of  acceptable 
values,  the  effective  tip  radii  are  supported  by  TEM  studies,  and  the  distribution  estimates  are  not  too 
dissimilar  from  experimental  expectations. 

Finally,  and  importantly,  note  that  by  constraining  the  statistical  distribution  to  be  given  by  one 
parameter  (%  tips  emitting,  AS,  or  a),  odd  results  may  occur.  For  example,  if  a  half-Gaussian 
distribution  for  an  array  for  which  only  a  subset  of  otherwise  highly  uniform  tips  are  working  is 
assumed,  a  large  c  value  will  be  estimated  in  order  to  lower  Afn  to  the  experimentally  observed  value. 
This  will  be  compensated  for  by  estimating  a  smaller  unit  cell  Bg  than  actually  characterizes  the 
subset.  In  the  absence  of  detailed  knowledge  of  the  distribution  of  tips,  such  problems  are 
unavoidable.  However,  we  note  that  for  the  experimental  data  considered,  the  width  of  the 
distribution  of  B’s  for  the  various  distributions  considered  are  similar,  arguing  against  such 
anomalous  situations. 

Further  improvement  and  vindication  of  the  modeling  effort  will  necessarily  include  information 
from  beam  diagnostics  (especially  with  regard  to  the  trajectory  analysis),  tip  characterization,  and 
array  scaling  behavior;  all  of  these  have  experimental  counterparts  within  the  NRL  program  in 
vacuum  microelectronics. 

Experimental  Results  vs  Boundary  Element  Model 

Figure  3  shows  the  discretization  of  the  conducting  surfaces  used  in  the  boundary  element  (BE) 
approach  for  the  Case  1  parameters  (the  Case  2  discretization  would  differ  only  in  scale).  Fine 
discretization  occurs  in  those  regions  where  the  field  variation  is  expected  to  be  large  (e.g.,  at  the  tip 
of  the  emitter,  near  gate  edges,  etc.).  The  parameters  used  in  the  BE  simulation  are  given  in  Table  2, 
and  correspond  to  the  labeling  in  Fig.  2.  The  results  of  the  simulation,  and  the  comparison  with 
experimental  data,  is  given  in  Table  2  in  the  section  “Statistical  Boundary  Element  Simulation.” 

Three  distributions  were  used  in  the  comparison  of  theory  with  the  cases  identified  in  Table  2: 

•  Delta  function  distribution:  tips  that  emitted  were  assumed  to  be  identical.  Consequently,  only 
one  value  of  Bfn  is  present  in  the  distribution.  The  number  of  tips  (or  %  tips  emitting)  was 
adjusted  to  achieve  correlation  between  theory  and  experiment  for  the  values  of  Afn  once 
was  identified. 
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Table  2  —  FEA  Unit  Cell  Simulation  Parameters 

Case  1  Case  2 

Symbol _ Designation _ (MIT) _ (SRI) _ Units 


Parameterization 


work  function 

2.0 

4.0 

eV 

th 

tip  height 

0.2 

1.3 

|im 

gate  radius 

0.075 

0.6 

|xm 

a 

tan"  ^  (base  radius/tip  height) 

20 

15 

degrees 

^8 

gate  thickness 

0.04 

0.3 

p^m 

dbg 

base  to  gate  distance 

0.195 

1.15 

pm 

dag 

gate  to  anode  distance 

0.5 

3.0 

pm 

N 

number  of  tips 

975 

1000 

- 

Va 

anode  voltage 

0.5 

50 

volts 

maximum  gate  voltage 

23 

78 

volts 

minimum  gate  voltage 

5 

39 

volts 

Statistical  boundary  element  simulation 

Afn(5) 

unit  cell 

318.7 

2.06 

MAA^2 

Bfn(6) 

unit  cell 

93.2 

423.8 

volts 

tip  radius 

30 

40 

A 

% 

percentage  tips  emitting 

65.85 

22.82 

% 

Afr(S) 

delta  distribution 

356.9 

5.6 

HAA^2 

Bfr(S) 

delta  distribution 

99.9 

467.2 

volts 

tip  radius 

33 

46 

A 

AB 

spread  in  B 

27 

300 

volts 

Afn  (L,) 

linear  distribution 

236.5 

2.27 

Bfn(L) 

linear  distribution 

100.5 

466.5 

volts 

C 

distribution  parameter 

25 

650 

volts 

Afn  (hG) 

half-Gaussian  distribution 

231.3 

1.17 

\xA/y^ 

BfnihG) 

half-Gaussian  distribution 

99.7 

466.0 

volts 

Experimental  parameterization  (least-squares  linear  fit)* 


^fn 

^fn 


235.0  1.28  MAA^2 

100.7  467  volts 
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Data  fit  used  a  weighted  least-squares  line  fit  on  the  experimental  points.  For  Case  I  (MIT):  9  V  <  <  15  V. 

For  Case  2  (SRI):  sample  53i+300-7Q  was  fitted  [9]. 
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•  Linear  distribution:  The  were  distributed  as  in  Eq.  (3).  The  value  of  AB  was  adjusted  to 
fix  Af„. 

•  Half-Gaussian:  The  were  distributed  as  in  Eq.  4,  where  B_  =  Bg.  The  value  of  ct  was 
adjusted  to  fix  Afn- 

The  consequences  of  these  three  distributions  with  regard  to  effective  tip  radius  are  indicated  in 
Table  2.  The  effective  tip  radius  was  assumed  to  be  the  same  for  both  the  unit  cell  and  the  linear  and 
half-Gaussian  distributions.  These  three  widely  disparate  assumptions  appear  to  result  in  comparable 
estimates  of  the  tip  radius,  and  in  the  similarity  of  the  AB  and  a  estimates  (the  full  width  at  half 
maximum  of  the  half-Gaussian  is  (V2)  a  ln(2),  which  is  comparable  to  AS).  Note  that  a  larger  gate 
radius  for  Case  2  engenders  a  larger  spread  in  the  S/„  values,  as  intuitively  expected.  Note  also  that 
the  similarity  of  effective  tip  radius  for  both  Case  1  and  Case  2  is  suportive  of  the  notion  of  surface 
nonuniformity.  These  estimates  are  comparible  to  the  findings  of  Levine  [8]  and  Cade  et  al.  [34]. 

From  the  distribution  in  Bf„'s  ,  we  can  infer  the  distribution  in  effective  tip  radius.  From  the 
tabulation  of  slopes  and  intercepts  in  the  semi-numerical  model  as  a  function  of  tip  radius,  we  can 
obtain  a  numerically  empirical  relation  between  the  B/>,’s  and  the  tip  radii,  analogous  to  Eq.  (38). 
For  brevity,  consider  the  results  for  Case  1  only.  Figure  4  shows  the  distribution  of  tips  entailed  by 
the  assumed  distribution  in  B’s.  The  area  under  the  linear  and  half-Gaussian  curves  is  equal  to  975 
(the  number  of  tips).  Note:  the  Delta  distribution  has  been  scaled  by  1/12,  and  there  is  a 
preponderance  of  tips  with  smaller  radii.  By  using  the  fitted  linear  relations  B/„  =  3 1 .23  -i-  2.078  a^ 
and  Afn  =-90.66  +  13.54  Oj,  in  volts,  and  |xA/v2,  respectively,  we  can  calculate  the  total  current 
contribution  from  each  group  of  tips  characterized  by  a  given  tip  radius,  as  shown  in  Fig.  5.  Notice 
that  the  total  current,  given  by  the  area  under  each  curve,  is  approximately  the  same.  Again,  compare 
the  findings  of  Cade  et  al.,  for  silicon  emitters  [34]. 

The  work  function  used  in  Figs.  4  and  5  was  for  a  cesiated  tip,  and  was  therefore  taken  to  be  2.0 
eV.  If  we  assume  that  the  geometry  and  distribution  of  tips  remains  identical,  but  we  now  consider 
the  tips  to  be  pure  molybdenum  with  a  work  function  of  4.35  eV,  we  can  predict  the  performance  of 
the  resulting  FEA  (Fig.  6).  On  an  FN  plot,  the  I{V)  characteristics  for  the  uncesiated  case  are,  to  a 
very  good  approximation,  linear,  and  described  by  Afn  =  22.90  |iA/V2  and  =  314  V. 

The  low  work  function  in  Case  1  raises  an  interesting  issue;  Is  approximating  J(F)  in  Eqs.  (27) 
and  (33)  by  Jfn(F)  valid?  Figure  7  compares  Jfn(F)  to  the  exact  (Airy  function  solution)  J(F)  [33]. 
For  typical  molybdenum  parameters  (and  at  room  temperature),  the  Fowler-Nordheim  estimate  is 
typically  too  low  by  20%  or  so  (the  discrepancy  becomes  worse  at  low  fields,  hence  voltages,  as  a 
result  of  Jfn(P)  being  a  zero-temperature  approximation).  For  the  cesiated  case,  Jfn(P)  significantly 
overestimates  the  actual  current  at  high  fields.  This  is  because  the  linearized  WKB  expression 
estimates  that  the  transmission  coefficient  is  greater  than  unity  in  the  region  where  tunneling  (and 
transmission  over  the  barrier)  is  important.  Using  the  correct  J(F)  in  the  expressions  for  I(V}  will 
result  in  a  bending  over  of  the  FN  curve  at  high  voltages  (low  V'^).  Consequently,  as  shown  in  Fig.  8, 
there  is  a  small  curvature  downward  at  low  (V'^)  when  using  the  exact  J(F).  Although  the  behavior  of 
the  experimental  points  in  this  regime  is  undoubtedly  primarily  due  to  space  charge  effects  [34], 
some  of  it  is  due  to  the  fact  that  nonlinearities  are  inherent  in  an  FN  plot  of  J(F}  in  this  regime. 
Finally,  note  that  the  exact  curve  is  somewhat  higher  than  that  using  JfnIP)'  as  a  result,  the  estimates 
of  the  spread  in  B  parameters  in  Table  2  are  smaller  than  would  otherwise  occur. 
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Tip  Radius  [Angstroms] 


Fig.  4  —  Distribution  of  emitters  according  to  tip  radius  for  Case  1  (MIT) 
parameters.  Total  number  of  emitters  =  975.  In  the  case  of  the  DeltaTunction 
distribution,  total  number  of  working  emitters  =  630 


Tip  Radius  [Angstroms] 

Fig.  5  —  Same  as  Fig.  4,  but  for  total  current  per  tip  radius  as  a  function  of  tip  radius 


Airy  FN^  ^  CuiTeilt  [|XA/tip] 
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Gate  Voltage  [Volts] 


—  Comparison  of  theoretical  cesiated  vs  uncesiated  FEAs  for  Case  1  (MIT) 
parameters  in  the  comparison  of  current  with  gate  voltage 


Fig.  7  —  Comparison  of  Fowler-Nordheim  approximation  to  J(F)  with  exact  solution 
for  Case  I  (MIT)  parameters  for  cesiated  and  uncesiated  work  functions 
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Fig.  8  —  Fowler-Nordheim  plot  of  experimental  data  and  boundary  element 
simulation.  Solid  line  is  for  the  calculation  using  the  exact  solution  of  J(F), 
dotted  line  is  for  the  calculation  using  the  FN  approximation  to  J(F).  Note  the 
curvature  at  high  gate  voltages  for  the  former 

Validity  of  the  Saturn  and  Semi-Numerical  Model  Approximations  as  Applied  to  the  Unit  Cell 

The  assumption  that  only  the  portion  of  the  tip  approximated  by  a  sphere  contributes  to  the  total 
current,  and  that  emission  along  the  shank  of  the  tip  is  negligible,  is  implicit  in  Eq.  (27).  In  fact,  the 
upper  limit  of  the  integral  in  Eq.  (26)  corresponds  to  an  angle  of  90°.  It  remains  to  be  demonstrated 
that  this  approximation  of  the  Saturn  model,  and  hence,  the  semi-numerical  model,  is  valid. 
Furthermore,  the  Saturn  model  itself  cannot  be  used  for  field-on-axis  estimates,  as  the  emitter  tip  is 
close  to  the  center  of  the  gate  plane  (see  Eq.  (25)  and  discussions  therein).  Similarly,  the  Saturn 
estimates  of  A  are  very  small,  whereas  for  the  BE  calculations,  A  =  0.32  for  a  range  of  geometries. 
Nevertheless,  the  semi-numerical  model  works  quite  well.  Figures  9  and  10  compare  the  full 
boundary  element  (BE)  calculation  for  the  parameters  listed  in  Table  2  with  the  semi-numerical  (SN) 
method,  in  which  the  linear  fits  shown  in  Eqs.  (36)  and  (37)  are  used  in  the  Saturn  model  equation 
for  current,  Eq.  (27).  As  can  be  seen,  the  estimates  are,  in  fact,  quite  close.  We  now  explore  the 
reasons  for  this.  Note  that  all  parameters  not  explicitly  mentioned  are  assumed  to  be  as  given  in 
Table  2  for  the  two  cases  considered. 

The  semi-numerical  method  implicitly  assumes  that  the  field  variation  along  the  surface  of  the 
emitter  is  governed  by  Eq.  (17),  albeit  that  the  and  A  parameter  are  given  by  the  boundary 
element  calculations,  and  are,  in  general,  gate  voltage-dependent.  F(6)  is  extracted  from  the 
boundary  element  data  by  defining  tanf0)  =  (Zs/Ps)<  where  6  is  the  polar  angle  of  the  equivalent 
sphere,  and  (p^,  Zs)  are  the  coordinates  along  the  surface  of  the  actual  emitter.  Figures  11  and  12 
show  the  F(d)  and  Ffp^)  comparisons  for  MIT-LL  and  SRI  geometries,  respectively,  as  given  in  Table 
2.  It  is  seen  that  for  "large"  angles,  the  deviation  between  the  BE  method  and  the  SN  method  grows. 
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Gate  Voltage  [Volts] 


Fig.  10  —  Comparison  of  the  boundary  element  and  the  semi-numerical  methods  (Case  2  parameters) 
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Fig.  1 1(a)  —  Field  along  the  surface  of  the  emitter  as  a  function  of  cylindrical 
coordinate  p  (Case  I  parameters  (MIT)) 


Fig.  1 1(b)  —  Field  along  the  surface  of  the  emitter  as  a  function  of  polar  angle  0 
(Case  1  parameters  (MIT)) 
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—  Field  along  the  surface  of  the  emitter  as  a  function  of  cylindrical  coordinate  p 
(Case  2  parameters  (SRI)) 


Fig.  12(b)  —  Field  along  the  surface  of  the  emitter  as  a  function  of  polar  angle  0 
(Case  2  parameters  (SRI)) 
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However,  the  current  density  falls  rapidly  with  field,  as  shown  in  Figs.  13  and  14,  and  hence,  the  SN 
approximation  is  good  despite  the  large  angle  behavior  (in  the  latter,  the  field  on  axis  tends  to  be 
mildly  noisy;  this  noise  can  be  diminished  with  proper  discretization,  and  is  in  progress).  The  field 
along  the  shank  of  the  emitter,  while  being  non-zero,  gives  rise  to  a  negligible  current,  and  therefore 
need  not  be  considered. 


Figures  15  and  16  show  the  behavior  of  the  5y„  and  inverse  field  enhancement  factor  (I/p)  as  a 
function  of  tip  and  gate  radius  for  the  two  cases  in  Table  2;  as  expected  from  Eq.  (38),  they  are 
predominantly  linear  in  tip  radius.  For  small  gate  radii,  the  linearity  in  In(ag)  is  approximate  and 
degrades  as  the  gate  radius  grows,  but  for  small  gate  radii  (where  the  Saturn  Model  is  presumed  to  be 
a  better  model),  the  linearity  predicted  by  Eq.  (38)  is  roughly  correct.  In  Figs  15(b)  and  16(b),  both 
(l/p)  and  Bfii  have  the  same  qualitative  dependence  on  gate  radius.  This  indicates  that  Bfn  is  linear  to 
a  good  approximation  in  (I/p),  as  can  be  shown  explicitly,  and  as  expected  from  Eq.  (38). 

Capacitance  of  FEA  Unit  Cell 


Figure  17  shows  the  behavior  of  the  ratio  of  Ctg  with  the  parallel  plate  approximation  given  by 


TTaj 


(39) 


based  on  the  dependencies  indicated  in  Eq.  (38).  The  x-axes  in  Fig.  17(a)  differ  from  those  in  Fig. 
17(b).  In  the  former,  the  tip  radius  is  varying,  while  in  the  latter,  it  is  not.  We  see  that  the  linearity 
prediced  by  Eq.  (38)  is  in  fact  maintained  in  the  boundary  element  calculations,  although  for  varying 
gate  radius,  the  linearity  is  apparent  when  the  gate  radius  is  small.  Recall  that  the  Saturn  model,  and 
hence,  Eq.  (38),  relies  on  the  assumption  that  the  gate  can  be  replaced  by  a  ring  of  charge,  and  that 
this  ring  lies  close  to  the  edge  of  the  gate.  When  the  gate  radius  becomes  large,  the  remainder  of  the 
gate,  and  the  charges  associated  with  it,  cannot  be  as  easily  dismissed  as  we  have  done. 

Note  that  what  we  are  examining  in  Fig.  17  is  the  tip-to-gate  capacitance  C,g,  not  the  capacitance 
of  the  FEA  unit  cell  overall.  Cfg  can  be  expected  to  be  dominated  by  the  parallel-plate-like 
capacitance  between  the  gate  and  the  base  if  the  tip-to-tip  separation  is  large.  The  charge  residing  on 
the  base  is  large;  in  fact,  the  percentage  of  the  charge  residing  on  the  tip  alone  of  the  total  tip+base 
charge  is  3.12%  for  Case  1  (MIT)  and  5.7%  for  Case  2  (SRI).  Consequently,  the  total  capacitance 
between  gate  and  (tip+base)  dominates 

Field  Distribution  of  FEA  Unit  Cell  for  Saturn  and  Boundary  Element  Models 

In  the  Saturn  model,  the  sphere  representing  the  tip  is  charged.  As  a  result,  far  from  the  sphere, 
the  field  due  to  that  charge  should  fall  off  in  a  1/r^  manner.  The  actual  charge  distribution  is 
responsible  for  higher  order  (multipole)  terms,  but  these  can  be  neglected.  In  the  boundary  element 
model,  two  effects  come  into  play:  (a)  the  charge  density  is  highest  on  the  emitter  tip,  and  (b)  the 
field  due  to  a  charged  ring  along  the  axis  of  the  ring  goes  as  the  inverse  square  of  the  distance  from 
the  observation  point  to  any  point  on  the  ring.  From  these  two  observations,  one  can  surmise  that 
near  the  emitter  tip,  the  field  fall-off  is  roughly  as  predicted  by  the  Saturn  model. 
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Fig.  13(a)  —  Current  density  along  the  surface  of  the  emitter  as  function  of  cylindrical 
coordinate  p  (Case  1  parameters  (MIT)) 
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Fig.  13(b)  —  Current  density  along  the  surface  of  the  emitter  as  function  of  polar  angle  0 

(Case  1  parameters  (MIT)) 
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Fig.  16(a)  —  Field  enhancement  factor  and  Fowler-Nordheim  B  value  as  function  of  tip  radius 

(Case  2  parameters) 
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Fig.  16(b)  —  Field  enhancement  factor  and  Fowler-Nordheim  B  value  as  function  of  gate  radius 

(Case  2  parameters) 
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Fig.  17(a)  —  Scaled  capacitance  for  Cases  1  and  2  parameters  (a)  tip  radius. 


Fig.  17(b)  —  Scaled  capacitance  for  Cases  1  and  2  parameters  (b)  field  enhancement. 
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An  expedient  method  for  calculating  both  the  potential  and  the  field  at  any  point  in  space  not  on 
a  conducting  surface  is  afforded  by  the  ring  representation  of  the  charged  surfaces,  as  in  Eq.  (34). 
On  axis  (p  =  0),  the  expressions  are  particularly  simple: 


mz)= 


Qi 


(40) 


Z-Zi 


Figure  18  shows  the  field  fall-off  as  calculated  by  Eq,  (40)  for  the  two  cases  considered  in  Table 
2.  Also  shown  is  the  Saturn  approximation,  which  can  be  roughly  recast  as 

Fsaturniz)  =  F,[^]^  ,  (41) 

where  Fq  is  the  field  at  the  tip,  and  as  is  the  radius  of  the  equivalent  sphere.  As  seen  in  Fig.  18,  Eq. 
(41)  accounts  for  most  of  the  behavior  near  the  emitter  tip. 

Incorporation  of  Semi-Numerical  Model  into  Larger  Simulations 

The  length  scale  characteristics  associated  with  the  unit  cell  FEA  analysis  preclude  the  boundary 
element  method,  as  presented  here,  from  being  used  for  an  entire  lOA  simulation,  although  it  is 
ideally  suited  for  modeling  conditions  near  the  tip.  Conversely,  techniques  using  a  finite-element 
analysis  [21-23]  are  ideally  suited  for  lOA  simulation,  whereas  the  discretization  required  for  atomic- 
dimension  tips  is  problematic.  A  method  of  combining  the  methods  is  indicated  by  observing  that 
the  boundary  element  method  requires  only  a  virtual  anode  to  model  conditions  near  the  emitter  tip. 
Consequently,  a  method  for  combining  boundary  element  with  finite-element  is  to  use  the  latter  to 
calculate  or  estimate  the  shape  and  potential  of  the  virtual  anode  and  then  use  the  former  to 
determine  emission  conditions.  The  field  conditions  in  the  immediate  vicinity  of  the  emitter  are  so 
extreme  that  space  charge  may  not  initially  hinder  the  implementation  of  this  algorithm,  although 
past  the  virtual  anode,  where  the  field  falls  off  considerably  and  the  characteristic  length  scales  are 
orders  of  magnitude  larger,  this  will  no  longer  in  general  be  true.  However,  in  the  post-virtual-anode 
region,  the  finite-element  method  is  ideally  suited  to  deal  with  space  charge  effects.  In  practice,  a  few 
iterations  between  a  finite-element  simulation  and  the  boundary  element  calculation  may  be  required 
to  correctly  include  the  effects  of  space  charge. 

CONCLUSIONS 

We  have  endeavored  to  present  a  simplified  model  of  the  FEA  unit  cell  (the  Saturn  model),  to 
show  that  the  analytical  formulae  suggested  by  this  model  are  valid,  and  to  show  how  this  model  can 
be  augmented  (the  semi-numerical  model)  in  order  to  correlate  with  a  full  numerical  solution  based 
on  the  boundary  element  method  of  the  field  emitter  unit  cell.  At  present,  attention  has  been 
restricted  to  refractory  metal  emitters  characterized  in  the  literature,  but  the  method  does  not  preclude 
an  extension  to  semiconductor  emitters.  We  have  identified  where  the  assumptions  in  both  the 
analytic  and  numerical  treatments  may  be  in  error,  and  have  shown  how  to  compensate  for  this. 
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Fig.  18(a)  —  Field  on  axis  (Case  I  (MIT)  parameters) 
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Fig.  18(b)  —  Field  on  axis  (Case  2  (SRI)  parameters) 
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Agreement  between  the  qualitative  behavior  of  the  actual  quantities  of  interest  (i.e.,  current-gate 
voltage  relations,  capacitance,  Fowler-Nordheim  Bfn  factors,  beta  factors,  field  and  potential 
calculations  within  the  unit  cell,  and  so  on)  and  the  Saturn  model  vindicates  the  latter's  simple 
representation  of  the  FEA  and  its  utility  in  understanding  FEA  unit  cell  behavior.  The  semi- 
numerical  method  demonstrates  how  the  quantitative  agreement  can  be  restored  between  Saturn  and 
boundary  element  through  the  introduction  of  slope  and  intercepts  for  Ftip  and  X  that  are  geometry 
dependent.  This  suggests  the  possibility  of  a  library  table-lookup  algorithm  to  treat  general  FEAs 
without  requiring  a  boundary  element  solution  of  the  problem.  Finally,  we  have  shown  how  the 
statistical  distribution  of  emitters  can  be  examined  to  generalize  from  unit  cell  analysis  to  arrays  in 
which  tip  non-uniformity  may  occur. 
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Appendix  A 

Linearity  of  I(V)  on  a  Fowler-Nordheim  Plot 

Consider  two  limiting  cases  of  the  current-voltage  relations  characterized  by  Eq.  (5).  The  “half- 
Gaussian”  or  “linear”  distribution  results  in  a  general  <I(V)>  relationship  of  the  form: 

/„(V)  =  v2  +  nexp(A^-5,/F).  (Ala) 

where  n  is  some  (not  necessarily  zero)  integer.  Note  that  while  corresponds  to  Ag  corresponds 
to  ln(A^).  A  Fowler-Nordheim  plot  of  /„(V),  in  which  ln(/„/V2)  corresponds  to  the  y-axis  and  1/V  to 
the  x-axis,  will  be  identically  linear  for  n  =  0,  with  a  slope  of  i-Bg)  and  a  y-intercept  of  Ag,  where  Ag 
and  Bg  are  constants.  The  “full-Gaussian”  distribution  results  in  a  general  <I(V)>  relationship  of  the 
form: 


W  =  v2exp 


(Alb) 


For  non-zero  n  and  cr,  the  question  arises  as  to  how  linear  ln(/„/V2)  or  ln(/o/v2)  will  appear. 


Polynomial  Least-Squares  Fitting  of  General  Functions 


A  function  y(x),  continuous  in  a  region  Xo  -  d  <x  <Xo  +  S,  is  expanded  in  terms  of  Legendre 
polynomials.  Let  ym(x)  be  that  polynomial  of  degree  m,  which  minimizes  the  L2  norm  of  PmM  - 
y(x),  where  Pm(x)  is  all  polynomials  of  degree  m  or  less  [35].  We  then  have 


yOT(jc)=  S' 
^  =  0 


where  z  =  (x-Xo)/8.  There  are  several  advantages  to  this  approach.  First,  the  Legendre  polynomials 
satisfy  -1  ^k(z)  <  1,  so  that  from  a  numerical  standpoint,  Eq.  (A2)  is  better  behaved  than  a  power 
series  expansion.  Second,  Pk(z)  satisfies  the  recursion  relation 

(k+l)Pk+i(z)  =  (2k+l)  z  Pkiz)  -  k  Pk.i(z)  (A3) 

which,  given  Po(z)  =  I  and  Pi(z)  =  z,  allows  the  higher  order  Legendre  polynomials  to  be  easily 
evaluated  for  a  given  z. 

The  coefficient  of  linear  correlation  R  gives  a  measure  of  how  well  y(x)  is  approximated  by 
ymM-  Let  us  define: 
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Then  R  =  Oxy  /  '^(<Jxx  <^yy)- 

Table  A1  -  Fowler-Nordheim  Linear  Fit  Parameters 


Parameter 

MIT-LL 

SRI 

Max  Vgate 

23  V 

78  V 

Min  Vgate 

5  V 

39  V 

Xo 

0.1217 

0.01923 

5 

0.0783 

0.00641 

Unit  Cell  Af^ 

318.7  ^AA^2 

2.063  pAA^2 

Unit  Cell  Bf^ 

93.18  V 

423.8  V 

Least-Squares  Fitting  of  Linear  or  half-Gaussian  y(x) 

For  the  particular  choice  y(x)  =  -  n  lr\(x/Xo)  +  A  -  J?  j:,  the  Ck  can  be  evaluated  analytically.  The 
first  few  have  been  found  to  be 


where 


Cq  =  A  - 
C,  =-B5+ 


o"  25 
3  n 


Co  = 


5  n 


12^3 


-x„S-i-[2-  f/)(5| 

Ixl-S^^-lx^S 
(52J5  +  2  2 <5^1  <5 


(A5) 


S  =  ln 


f 

s 

2^ 

;  U=ln 

1  - 

Xo 

1 

(A6) 


The  linearity  of  y(x)  is  indicated  by  how  close  is  to  unity,  or  by  the  smallness  of  C2  by 
comparison  to  Cq  and  C/.  Finally,  for  a  linear  fit  (n  =  /),  in  which  y](x)  =  A  -  5  a:,  the  best  A  and  B 
are  found  to  be 
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(A7) 


Consider  the  two  experimentally  motivated  cases  discussed  earlier,  denoted  MIT-LL  [Al,  A2]  and 
SRI  [A3].  Table  Al  shows  the  relevant  parameters,  based  on  a  least-squares  fit  of  the  experimental 
data  and  its  relation  to  the  unit  cell  modeling  previously  described.  Tables  2  and  A3  show  the  results 
of  the  best  linear  fit  calculations.  As  indicated  by  both  the  value  of  C2  and  R,  even  for  n  large,  y„(x) 
is  still,  to  a  good  approximation,  linear.  A  positive  value  of  C2  indicates  that  the  curve  is  concave-up. 

Table  A2  -  Fowler-Nordheim  Linearity  of  /^(VO 


n 

Afy,  [nM"'2} 

Bfn  [Volts] 

Co 

Cl 

C2 

R 

Case  1 

;  MIT-LL  Unit  ^ 

Cell 

0 

318.7 

93.2 

-19.40 

-7.293 

0.000 

1.0000 

1 

1039 

102.2 

-19.32 

-8.001 

0.1698 

0.9999 

2 

3386 

111.3 

-19.24 

-8.709 

0.3397 

0.9995 

3 

1.104e+04 

120.3 

-19.16 

-9.418 

0.5095 

0.9991 

4 

3.597e+04 

129.4 

-19.08 

-10.13 

0.6794 

0.9986 

5 

1.172e+05 

138.4 

-19.00 

-10.83 

0.8492 

0.9980 

Case  2:  SRI  Unit  Cell 

0 

2.063 

423.8 

-21.24 

-2.717 

0.000 

1.0000 

1 

5.850 

477.0 

-21.22 

-3.058 

0.03892 

1.0000 

2 

16.59 

530.3 

-21.20 

-3.399 

0.07784 

0.9998 

3 

47.06 

583.5 

-21.18 

-3.740 

0.1168 

0.9997 

4 

133.5 

636.7 

-21.17 

-4.081 

0.1557 

0.9996 

5 

378.6 

689.9 

-21.15 

-4.422 

0.1946 

0.9994 

Least-Squares  Fitting  of  full-Gaussian  y(x) 

For  the  particular  choice  y(x)  =  Ao~  Bg  x  +  D  x^,  where  D  =  6^/4,  the  Ck  can  also  be  evaluated 
analytically:  for  k>  2,  the  Ck  vanish.  The  remainder  are  given  by 

CQ=Ao-B^o  +  \D{5^  +  ^x^^  (A8) 

C^={-B„  +  2Dx„)5 
C2  =  \Dd^  . 


In  contrast  to  Eq.  (A6),  non-zero  o  will  cause  to  decrease;  the  equation  for  y{x)  is  a  concave-up 
parabola.  For  the  same  parameters  in  Table  Al,  we  can  construct  in  Table  A3  an  analog  to  Table  A2. 
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Table  A3  -  Fowler-Nordheim  Linearity  of  I(^V) 


(y[V]  Afn  [4AA^^2]  Bf„  fVoltsl 

Co 

C} 

Cl  R  1 

Case  I:  MIT-LL  Unit  Cell 


318.7 

93.18 

-19.39 

-7.292 

.0000 

294.2 

91.66 

-19.29 

-7.173 

.0255 

231.5 

87.09 

-18.97 

-6.816 

.1021 

155.3 

79.48 

-18.45 

-6.221 

.2297 

88.80 

68.83 

-17.71 

-5.387 

.4083 

43.27 

55.14 

-16.76 

-4.315 

.6380 

.0 

5.0 

10.0 

15.0 

20.0 

25.0 


Case!;  SRI  Unit  Cell 


.0 

30.0 

60.0 

90.0 

120.0 

150.0 


2.063 

423.8 

-21.24 

-2.717 

.0000 

1.904 

415.1 

-21.16 

-2.661 

0.0061 

1.497 

389.2 

-20.90 

-2.495 

0.0247 

1.003 

345.9 

-20.46 

-2.217 

0.0555 

.5724 

285.3 

-19.86 

-1.829 

0.0986 

.2783 

207.5 

-19.08 

-1.330 

0.1541 

1.0000 

1.0000 

0.9999 

0.9996 

0.9983 

0.9935 


1.0000 

1.0000 

1.0000 

0.9998 

0.9991 

0.9960 


In  this  case,  for  a  given  Ag  and  Bo,  the  effects  of  a  wider  statistical  distribution  of  B’s,  as 
represented  by  a,  results  in  a  smaller  best  fit  Bf„  value.  This  is  in  constrast  to  the  linear  or  half- 
Gaussian  results.  However,  note  that  in  the  former  case.  Bo  represents  the  smallest  B  value  in  the 
distribution,  whereas  in  the  latter,  it  represents  the  mean  value. 
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Appendix  B 

Chemical  Potential  and  Temperature-Dependent  Fowler-Nordheim  Equation 


The  chemical  potential  and  temperature-dependent  Fowler-Nordheim  equation  for  metallic 
parameters  is  given  by 


1 


y/p) 

n(cyh7t/P) 


CLfn  — 


^  167t^h(t>r(y)^ 


bfr  =  ^\p2^v(y) 


<^yh=:^N/2m<t»  /(y)  , 


(Bl) 


where  y'^  =  4QF/<t)^.  Note  that  and  as  defined  here,  are  not  independent  of  the  field  F.  The 
functions  v(y)  and  r(y)  can  be  evaluated  via  integral  definitions.  Polynomials  appropriate  for 
numerical  work  are  [Bl]: 


t{y) 

v(y) 

ziy) 


n 


|3  { [(4298z  +  465)z  +  40]  z  +  4}  z  -  8 J  z  +  4 


^  1 
Z-fl 

=  -^  71  ({  -  [3(665z  +  94)z  +  52]z  -  20}z  +  4)z 
■  16  • 


(B2) 


It  is  common  in  the  literature  to  take  the  zero  temperature  {T  — >  0)  and  infinite  chemical  potential 
(^i  oo)  limit  of  Eq.  (Bl)  (which  amounts  to  neglecting  the  terms  in  parentheses),  and  approximate 
f(y)  by  a  constant  and  v(y)  as  a  linear  function  in  y^.  By  Taylor-expanding  v(y)  about  y„  and 
insisting  that  the  coefficient  of  y^  be  unity,  it  can  be  shown  that 

t(y)  =  1.05657 
v(y)  =  0.93685  -y^ 

Equation  (B3),  however,  is  not  the  least-squares  approximation,  which  can  be  shown  to  be: 


t(y)  =  1.06667  (B4) 

v(y)  =  0.94199  -  0.96969y  2 

In  spite  of  the  fact  that  Eq.  (B4)  is  closer  to  the  actual  values  of  v(y)  and  t(y),  the  form  of  Eq.  (B3)  is 
often  used  because  of  the  simplicity  of  the  y^  coefficient. 


Appendix  C 

Field  Enhancement  Effects  Due  To  Surface  Curvature 


Boss  on  Plane 

Consider  an  infinite  planar  conductor  with  a  hemispherical  boss  of  radius  a  on  the  surface.  Far 
from  the  conductor,  let  the  field  F  in  the  z  direction  be  constant.  Using  a  Legendre  polynomial 
expansion  method  [Cl],  the  potential  everywhere  is 


V(r,0)  =V^-Frcos(0)[l-(f)^ 


(Cl) 


where  z  =  r  cos(0).  The  field  along  the  surface  of  the  boss  (r  =  a)  is  then 

F{a,d)  =  -  VV(r,0)  =  3  Fcos  (0).  (C2) 

The  field  at  the  tip  of  a  protrusion  is  therefore  enhanced  by  a  factor  of  3  compared  to  the  field  far 
away  from  the  boss. 

Two-Dimensional  Hyperbolic  Wedge  (Diode) 


For  generalized  coordinates  {qhq2y(l3\  the  gradient  and  the  Laplacian  are  given  by  [C2] 


hj  =  djx  +  djy  +  djz 

V2  =  (/i  JA2/13)  ‘  ^  ^(1  +  £ijk)£ijk 

where  is  the  Levi-Civita  symbol,  and  {qi,q2,(l3)  =  (x,y,z)  in  Cartesian  coordinates, 
dimensional  hyperbolic  coordinates  (a,p,y)  we  have: 


(C3) 


In  2- 


x  =  a^^  sinh  (a)  sin  {p) 

y  =  aij  cosh  (a)  cos  (p)  (C4) 

z  =  y. 


where  Oh  is  the  distance  from  the  origin  to  the  focus  of  the  hyperbola  defined  by  constant  p.  Solving 
V^y/=  0  under  the  boundary  conditions  y/(P=Po)  =  0  (“wedge”)  and  \if(P=7C/2)  =  Vq  (“anode”)  is  a 
straightforward  problem  and  results  in 
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_  Vq 

^h{i  - Po]  V'sinh^(a)  +  sin2(y?)' 


(C5) 


Define  the  radius  of  curvature  of  the  wedge  as  y  l;c=o  =  1/a^-  and  the  distance  from  the  wedge 
to  the  anode  as  such  that  tan^  /  Zo-  In  the  limit  that  Po  is  small,  the  field  at  the  tip  goes  as 

F'tip  ^2Vo/  ivJ(as  Zo). 


Three-Dimensional  Hyperbolic  Tip  (Diode) 


Consider  a  hyperbolic  tip  with  rotational  symmetry,  for  which 

x  =  af^  sinh  (a)  sin  {p)  cos  (f) 

y  =  af,sinh(a)sm(p)sm(y)  (C6) 

z  =  afj  cosh  (or)  cos  (p) . 

Invoking  Eq.  (C3),  letting  y/=  X(a)  Y(^  Z(y),  and  exploiting  the  independence  of  i//on  ygives 


3„^sinh  (a)  -  n(n  +  1)  sinh  (a) 

3yg(sin  (P)  d^j  +  n(n  +  1)  sin(p) 


X(a)  =  0 

y(P)=o 


(C7) 


for  which 


X„(ar)  =  a  P„(cos  (la))  +  a'  e„(cos  (/or)) 

Y^ip)  =  b  P„(cos  (p)  +  b'  e„(cos  ip) , 

where  P „(x)  and  Q„(x)  are  Legendre  polynomials  of  the  first  and  second  kind,  the  a's  and  b's  are 
constants.  Invoking  the  boundary  conditions,  analogous  to  the  wedge,  that  the  tip  is  at  zero  potential 
and  the  anode  at  we  find  that  n  =  0  and 


J  go(cos  {p) 


l2o(cOS  {p^)] 


Fia,p  =  (-p- 


j2o(cos  (^p)j  a/j  sin  (P  ^sinh^for)  +  sin^(y8) 
In  the  limit  of  small  the  field  at  the  tip  goes  as 


(C9) 


(CIO) 


Note  that  this  is  a  slightly  weaker  dependence  on  the  tip  radius  than  in  Eq.  (38d)  for  the  triode,  as  the 
In-term  in  Eq.  (CIO)  contains  a  factor  of  a^.  Invoking  the  polar  coordinate-based  parametric 
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representation  introduced  in  Eq.  (17).  we  find  that  A=  l/cosHPo)  =  1-  Equation  (27)  may  now  be 
used  to  estimate  the  current  from  a  hyperbolic  diode. 

REFERENCE 

Cl.  G.  Arfken,  Mathematical  Methods  for  Physicists,  3rd  edition  (Academic  Press,  Orlando,  Florida, 
1985). 


